Invariant grids for reaction kinetics 



Alexander N. Gorban 

Department of Materials, Institute of Polymer Physics 
Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland 
Institute of Computational Modeling RAS, Krasnoyarsk 660036, Russia 

Iliya V. Karlin 

Department of Materials, Institute of Polymer Physics 
Swiss Federal Institute of Technology, CH-8092 Zurich, Switzerland 

Andrei Yu. Zinovyev 

Institut des Hautes Etudes Scientifiques, 
Le Bois-Marie, 35, route de Chartres, F-91440, Bures-sur-Yvette, France 



Abstract 

In this paper, we review the construction of low-dimensional manifolds of reduced 
description for equations of chemical kinetics from the standpoint of the method 
of invariant manifold (MIM). MIM is based on a formulation of the condition of 
invariance as an equation, and its solution by Newton iterations. A grid-based ver- 
sion of MIM is developed. Generalizations to open systems are suggested. The set 
of methods covered makes it possible to effectively reduce description in chemical 
kinetics. 

The most essential new element of this paper is the systematic consideration of a 
discrete analogue of the slow (stable) positively invariant manifolds for dissipative 
systems, invariant grids. We describe the Newton method and the relaxation method 
for the invariant grids construction. The problem of the grid correction is fully 
decomposed into the problems of the grid's nodes correction. The edges between the 
nodes appears only in the calculation of the tangent spaces. This fact determines 
high computational efficiency of the invariant grids method. 
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1 Introduction 



In this paper, we present a general method of constructing the reduced description for 
dissipative systems of reaction kinetics and a new method of invariant grids. Our approach 
is based on the method of invariant manifold which was developed in end of 1980th - be- 
ginning of 1990th [24,25,26]. Its realization for a generic dissipative systems was discussed 
in [28,36]. This method was applied to a set of problems of classical kinetic theory based 
on the Boltzmann kinetic equation [28,49,51]. The method of invariant manifold was suc- 
cessfully applied to a derivation of reduced description for kinetic equations of polymeric 
solutions [79]. It was also been tested on systems of chemical kinetics [35,32]. In order to 
construct manifolds of a relatively low dimension, grid-based representations of manifolds 
become a relevant option. The idea of invariant grids was suggested recently in [32]. 

The goal of nonequilibrium statistical physics is the understanding of how a system with 
many degrees of freedom acquires a description with a few degrees of freedom. This 
should lead to reliable methods of extracting the macroscopic description from a detailed 
microscopic description. 

Meanwhile this general problem is still far from the final solution, it is reasonable to 
study simplified models, where, on the one hand, a detailed description is accessible to 
numerics, on the other hand, analytical methods designed to the solution of problems in 
real systems can be tested. 

In this paper we address the well known class of finite-dimensional systems known from 
the theory of reaction kinetics. These are equations governing a complex relaxation in 
perfectly stirred closed chemically active mixtures. Dissipative properties of such systems 
are characterized with a global convex Lyapunov function G (thermodynamic potential) 
which implements the second law of thermodynamics: As the time t tends to infinity, 
the system reaches the unique equilibrium state while in the course of the transition the 
Lyapunov function decreases monotonically. 

While the limiting behavior of the dissipative systems just described is certainly very 
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simple, there are still interesting questions to be asked about. One of these questions is 
closely related to the above general problem of nonequilibrium statistical physics. Indeed, 
evidence of numerical integration of such systems often demonstrates that the relaxation 
has a certain geometrical structure in the phase space. Namely, typical individual tra- 
jectories tend to manifolds of lower dimension, and further proceed to the equilibrium 
essentially along these manifolds. Thus, such systems demonstrate a dimensional reduc- 
tion, and therefore establish a more macroscopic description after some time since the 
beginning of the relaxation. 

There are two intuitive ideas behind our approach, and we shall now discuss them in- 
formally. Objects to be considered below are manifolds (surfaces) Q in the phase space 
of the reaction kinetic system (the phase space is usually a convex polytope in a finite- 
dimensional real space). The 'ideal' picture of the reduced description we have in mind 
is as follows: A typical phase trajectory, c(t), where t is the time, and c is an element 
of the phase space, consists of two pronounced segments. The first segment connects the 
beginning of the trajectory, c(0), with a certain point, c(ti), on the manifold Q (rigorously 
speaking, we should think of c(ti) not on f2 but in a small neighborhood of fl but this 
is inessential for the ideal picture). The second segment belongs to fl, and connects the 
point c(t±) with the equilibrium c cq = c(oo), c cq e f2. Thus, the manifolds appearing in 
our ideal picture are "patterns" formed by the segments of individual trajectories, and 
the goal of the reduced description is to "filter out" this manifold. 

There are two important features behind this ideal picture. The first feature is the in- 
variance of the manifold Q: Once the individual trajectory has started on f2, it does not 
leaves Q anymore. The second feature is the projecting: The phase points outside ft will 
be projected onto ft. Furthermore, the dissipativity of the system provides an additional 
information about this ideal picture: Regardless of what happens on the manifold f2, the 
function G was decreasing along each individual trajectory before it reached f2. This ideal 
picture is the guide to extract slow invariant manifolds. 

One more point needs a clarification before going any further. Low dimensional invariant 
manifolds exist also for systems with a more complicated dynamic behavior, so why to 
study the invariant manifolds of slow motions for a particular class of purely dissipative 
systems? The answer is in the following: Most of the physically significant models in- 
clude non-dissipative components in a form of either a conservative dynamics, or in the 
form of external forcing or external fluxes. Example of the first kind is the free flight of 
particles on top of the dissipation-producing collisions in the Boltzmann equation. For 
the second type of example one can think of irreversible reactions among the suggested 
stoichiometric mechanism (inverse process are so unprobable that we discard them com- 
pletely thereby effectively "opening" the system to the remaining irreversible flux). For 
all such systems, the present method is applicable almost without special refinements, 
and bears the significance that invariant manifolds are constructed as a "deformation" of 
the relevant manifolds of slow motion of the purely dissipative dynamics. Example of this 
construction for open systems is presented below in section 11. Till then we focus on the 
purely dissipative case for the reason just clarified. 

The most essential new element of this paper is the systematic consideration of a dis- 
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crete analogue of the slow (stable) positively invariant manifolds for dissipative systems, 
invariant grids. These invariant grids were introduced in the [32]. Here we will describe 
the Newton method subject to incomplete linearization and the relaxation methods for 
the invariant grids. It is worth to mention, that the problem of the grid correction is 
fully decomposed into the problems of the grid's nodes correction. The edges between the 
nodes appears only in the calculation of the tangent spaces. This fact determines high 
computational efficiency of the invariant grids method. 

Due to the famous Lyapunov auxiliary theorem [60,54] we can construct analytical in- 
variant manifolds for kinetic equations with analytical right hand side. Moreover, the 
analycity can serve as a "selection rule" for selection the unique analytic positively invari- 
ant manifold from the infinite set of smooth positively invariant manifolds. The analycity 
gives a possibility to use the powerful technique of analytical continuation and Carleman's 
formulae [1,38,39]. It leads us to superresolution effects: A small grid may be sufficient to 
present an "large" analytical manifold immersed in the whole space. 

The paper is organized as follows. In the section 2, we review the reaction kinetics (section 
2.1), and discuss the main methods of model reduction in chemical kinetics (section 2.2). In 
particular, we present two general versions of extending partially equilibrium manifolds to 
a single relaxation time model in the whole phase space, and develop a thermodynamically 
consistent version of the intrinsic low-dimensional manifold (ILDM) approach. In the 
section 3 we review the method of invariant manifold in the way appropriate to this 
class of nonequilibrium systems. In the sections 4 and 5 we give some details on the 
two relatively independent parts of the method, the thermodynamic projector, and the 
iterations for solving the invariance equation. 

We also describe a general symmetric linearization procedure for the invariance equation, 
and discuss its relevance to the picture of decomposition of motions. In the section 6, these 
two procedures are combined into an unique algorithm. In the section 7, we demonstrate an 
illustrative example of analytic computations for a model catalytic reaction. In the section 
8 we introduce the relaxation method for solution the invariance equation. This relaxation 
method is an alternative to the Newton iteration method. In the section 9 we demonstrate 
how the thermodynamic projector is constructed without the a priori parameterization 
of the manifold 1 . This result is essentially used in the section 10 where we introduce 
a computationally effective grid-based method to construct invariant manifolds. It is the 
central section of the paper. We present the Newton method and the relaxation method for 
the grid construction. The Carleman formulas for analytical continuation a manifold from 
a grid are proposed. Two examples of kinetic equations are analyzed: a two-dimensional 
catalytic reaction (four species, two balances) and a four-dimensional oxidation reaction 
(six species, two balances). 

In the section 11 we describe an extension of the method of invariant manifold to open 



This thermodynamic projector is the unique operator which transforms the arbitrary vector 
field equipped with the given Lyapunov function into a vector field with the same Lyapunov 
function (and also this happens on any manifold which is not tangent to the level of the Lyapunov 
function). 
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systems. Finally, results are discussed in the section 12. 



2 Equations of chemical kinetics and their reduction 

2. 1 Outline of the dissipative reaction kinetics 

We begin with an outline of the reaction kinetics (for details see e. g. the book of [77]). Let 
us consider a closed system with n chemical species Ai, . . . , A n , participating in a complex 
reaction. The complex reaction is represented by the following stoichiometric mechanism: 

a,iAi + . . . + a sn A n ^ /3 al Ai + . . . + (3 sn A n , (1) 

where the index s — 1, . . . ,r enumerates the reaction steps, and where integers, a S i and 
/3 S i, are stoichiometric coefficients. For each reaction step s, we introduce n-component 
vectors a s and j3 s with components a S i and (3 si . Notation 7 S stands for the vector with 
integer components 7 S j = (5 S i — a S i (the stoichiometric vector). We adopt an abbreviated 
notation for the standard scalar product of the n-component vectors: 

n 

i=i 

The system is described by the n-component concentration vector c, where the component 
Cj > represents the concentration of the specie Aj. Conservation laws impose linear 
constraints on admissible vectors c (balances): 

(b i ,c)=B i , i = l,...,l, (2) 

where bi are fixed and linearly independent vectors, and Bi are given scalars. Let us denote 
as B the set of vectors which satisfy the conservation laws (2): 

B = {c\(b 1 ,c)=B 1 ,...,(b l ,c) = B l }. 

The phase space V of the system is the intersection of the cone of n-dimensional vectors 
with nonnegative components, with the set B, and dimV = d — n — I. In the sequel, we 
term a vector c e V the state of the system. In addition, we assume that each of the 
conservation laws is supported by each elementary reaction step, that is 

h s ,b l ) = 0, (3) 
for each pair of vectors 7 S and bi. 
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Reaction kinetic equations describe variations of the states in time. Given the stoichio- 
metric mechanism (1), the reaction kinetic equations read: 

c= J(c), J(c) = 5>,W.(c), (4) 

3 = 1 



where dot denotes the time derivative, and W s is the reaction rate function of the step s. 
In particular, the mass action law suggests the polynomial form of the reaction rates: 



w s = kt\{cT-K\{4% (5) 



where kf and k~ are the constants of the direct and of the inverse reactions rates of the 
sth reaction step. The phase space V is positive-invariant of the system (4): If c(0) G V, 
then c(t) E V for all the times t > 0. 

In the sequel, we assume that the kinetic equation (4) describes evolution towards the 
unique equilibrium state, c eq , in the interior of the phase space V. Furthermore, we assume 
that there exists a strictly convex function G(c) which decreases monotonically in time 
due to Eq. (4) 2 : 

G = (VG(c),J(c))<0, (6) 



Here VG is the vector of partial derivatives dG/dci, and the convexity assumes that the 
n x n matrices 

H c =\\d 2 G(c)/dc l dc J \\, (7) 



are positive definite for all c e V. In addition, we assume that the matrices (7) are 
invertible if c is taken in the interior of the phase space. 

The function G is the Lyapunov function of the system (4), and c eq is the point of global 
minimum of the function G in the phase space V. Otherwise stated, the manifold of 
equilibrium states c eq (i?i, . . . ,Bi) is the solution to the variational problem, 

G — > min for (bj, c) = i — 1, . . . , I. (8) 

2 With some abuse of language, we can term the functional — G the entropy, although it is a 
different functional for non-isolated systems. We recall that thermodynamic Lyapunov functions 
are well defined not just for isolated systems. Such functional are easily constructed also for 
systems which exchange energy and/or matter with a larger equilibrium system (with a thermo- 
stat, for example). In such a case, the thermodynamic Lyapunov function is constructed as the 
entropy of the minimal closed system containing the system under consideration [23] . In partic- 
ular, the free energy and the free enthalpy (the Gibbs and the Helmholz energies, respectively) 
can be constructed in this manner. They are are identical with the entropy of the minimal closed 
system containing the given system within the accuracy of multiplication with a factor which 
remains constant in time, and subtracting a constant. 



7 



For each fixed value of the conserved quantities B{, the solution is unique. In many cases, 
however, it is convenient to consider the whole equilibrium manifold, keeping the conserved 
quantities as parameters. 

For example, for perfect systems in a constant volume under a constant temperature, the 
Lyapunov function G reads: 

n 

G = E^NQ/cD-l]. (9) 
i=i 



It is important to stress that c cq in Eq. (9) is an arbitrary equilibrium of the system, under 
arbitrary values of the balances. In order to compute G(c), it is unnecessary to calculate 
the specific equilibrium c cq which corresponds to the initial state c. Moreover, for ideal 
systems, function G is constructed from the thermodynamic data of individual species, 
and, as the result of this construction, it turns out that it has the form of Eq. (9). Let us 
mention here the classical formula for the free energy F = RTVG: 

n 

f = VRTJ2a[(HaVQ *) - i) + *CH], (10) 
i=i 



where V is the volume of the system, T is the temperature, Vq j = NQ(2irh 2 /mikT) 3 / 2 
is the quantum volume of one mole of the specie Ai, Nq is the Avogadro number, is 
the mass of the molecule of Aj, R = kNo, and Fi nt j(T) is the free energy of the internal 
degrees of freedom per mole of Aj. 

Finally, we recall an important generalization of the mass action law (5), known as the 
Marcelin-De Donder kinetic function. This generalization was developed in [20] based on 
ideas of the thermodynamic theory of affinity [18]. We use the kinetic function suggested 
in [11]. Within this approach, the functions W s are constructed as follows: For a given 
strictly convex function G, and for a given stoichiometric mechanism (1), we define the 
gain (+) and the loss (— ) rates of the sth step, 

W, + = ^exp[(VG,a.)], W s ~ = <p; exp[(VG, /3J], (11) 



where ipf > are kinetic factors. The Marcelin-De Donder kinetic function reads: W s = 
— W~ , and the right hand side of the kinetic equation (4) becomes, 

J = jl fsWt exp[(VG, a.)] - y?; exp[(VG, (3 S )]}. (12) 

s=l 

For the Marcelin-De Donder reaction rate (11), the dissipation inequality (6) reads: 

G = £[(VG, p.) - (VG, a.)] Ue {VG ' a ° ] - <P7e^ G M < 0. (13) 

s=l > 
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The kinetic factors Lpf should satisfy certain conditions in order to make valid the dissi- 
pation inequality (13). A well known sufficient condition is the detail balance: 

V>t = <P7, (14) 

other sufficient conditions are discussed in detail elsewhere [77,23,46,47]. For the function 
G of the form (9), the Marcelin-De Donder equation casts into the more familiar mass 
action law form (5). 

2.2 The problem of reduced description in chemical kinetics 

What does it mean, "to reduce the description of a chemical system"? This means the 
following: 

(1) To shorten the list of species. This, in turn, can be achieved in two ways: 

(i) To eliminate inessential components from the list; 

(ii) To lump some of the species into integrated components. 

(2) To shorten the list of reactions. This also can be done in several ways: 

(i) To eliminate inessential reactions, those which do not significantly influence the 
reaction process; 

(ii) To assume that some of the reactions "have been already completed" , and that 
the equilibrium has been reached along their paths (this leads to dimensional reduc- 
tion because the rate constants of the "completed" reactions are not used thereafter, 
what one needs are equilibrium constants only). 

(3) To decompose the motions into fast and slow, into independent (almost-independent) 
and slaved etc. As the result of such a decomposition, the system admits a study "in 
parts". After that, results of this study are combined into a joint picture. There are 
several approaches which fall into this category: The famous method of the quasi- 
steady state (QSS), pioneered by Bodenstein and Semenov and explored in consider- 
able detail by many authors, in particular, in [10,14,68,21,66], and many others; the 
quasi-equilibrium approximation [62,23,74,21,46,47]; methods of sensitivity analysis 
[64,56]; methods based on the derivation of the so-called intrinsic low-dimensional 
manifolds (ILDM, as suggested in [61]). Our method of invariant manifold (MIM, 
[24,25,26,28,35,36]) also belongs to this kind of methods. 

Why to reduce description in the times of supercomputers? 

First, in order to gain understanding. In the process of reducing the description one is 
often able to extract the essential, and the mechanisms of the processes under study 
become more transparent. Second, if one is given the detailed description of the system, 
then one should be able also to solve the initial-value problem for this system. But what 
should one do in the case where the the system is representing just a point in a three- 
dimensional flow? The problem of reduction becomes particularly important for modeling 
the spatially distributed physical and chemical processes. Third, without reducing the 
kinetic model, it is impossible to construct this model. This statement seems paradoxal 
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only at the first glance: How come, the model is first simplified, and is constructed only 
after the simplification is done? However, in practice, the typical for a mathematician 
statement of the problem, (Let the system of differential equations be given, then ...) 
is rather rarely applicable in the chemical engineering science for detailed kinetics. Some 
reactions are known precisely, some other - only hypothetically. Some intermediate species 
are well studied, some others - not, it is not known much about them. Situation is even 
worse with the reaction rates. Quite on the contrary, the thermodynamic data (energies, 
enthalpies, entropies, chemical potentials etc) for sufficiently rarefied systems are quite 
reliable. Final identification of the model is always done on the basis of comparison with 
the experiment and with a help of fitting. For this purpose, it is extremely important to 
reduce the dimension of the system, and to reduce the number of tunable parameters. 
The normal logics of modeling for the purpose of chemical engineering science is the 
following: Exceedingly detailed but coarse with respect to parameters system — > reduction 
— > fitting — > reduced model with specified parameters (cycles are allowed in this scheme, 
with returns from fitting to more detailed models etc). A more radical viewpoint is also 
possible: In the chemical engineering science, detailed kinetics is impossible, useless, and 
it docs not exist. For a recently published discussion on this topic see [57,58]; [76]. 

Alas, with a mathematical statement of the problem related to reduction, we all have 
to begin with the usual: Let the system of differential equations be given .... Enormous 
difficulties related to the question of how well the original system is modeling the real 
kinetics remain out of focus of these studies. 

Our present work is devoted to studying reductions in a given system of kinetic equa- 
tions to invariant manifolds of slow motions. We begin with a brief discussion of existing 
approaches. 



2.3 Partial equilibrium approximations 



Quasi- equilibrium with respect to reactions is constructed as follows: From the list of 
reactions (1), one selects those which are assumed to equilibrate first. Let they be indexed 
with the numbers s±, . . . , s^. The quasi-equilibrium manifold is defined by the system of 
equations, 

W+ = W~, < = !,...,*. (15) 



This system of equations looks particularly elegant when written in terms of conjugated 
(dual) variables, /x = VG: 

(7..,/z) = 0, 1 = 1,...,*. (16) 



In terms of conjugated variables, the quasi-equilibrium manifold forms a linear subspace. 
This subspace, L- 1 , is the orthogonal completement to the linear envelope of vectors, 
£ = lm {7 Sl ,---,7 s J- 
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Quasi- equilibrium with respect to species is constructed practically in the same way but 
without selecting the subset of reactions. For a given set of species, A^ 17 . . . , ^4j fe , one 
assumes that they evolve fast to equilibrium, and remain there. Formally, this means 
that in the /c-dimensional subspace of the space of concentrations with the coordinates 
Cjj , . . . , Ci k , one constructs the subspace L which is defined by the balance equations, 
(bi,c) = 0. In terms of the conjugated variables, the quasi-equilibrium manifold, L- 1 , is 
defined by equations, 

/xGL x , (/i= (17) 

The same quasi-equilibrium manifold can be also defined with the help of fictitious reac- 
tions: Let g 1: . . . , g q be a basis in L. Then Eq. (17) may be rewritten as follows: 

(^,/x) = 0, i = l,...,q. (18) 

Illustration: Quasi-equilibrium with respect to reactions in hydrogen oxidation: Let us 
assume equilibrium with respect to dissociation reactions, H 2 2H, and, 2 ^ 20, in 
some subdomain of reaction conditions. This gives: 

k± c H 2 = ki C H , fcj c 2 = ^2 c Cr 

Quasi-equilibrium with respect to species: For the same reaction, let us assume equilibrium 
over H, O, OH, and H 2 2 , in a subdomain of reaction conditions. Subspace L is defined 
by balance constraints: 

ch + coh + 2ch 2 o 2 =0, c + c H + 2ch 2 o 2 = 0. 

Subspace L is two-dimensional. Its basis, {gi,g 2 } m the coordinates ch, Co, coh, an d 
ch 2 o 2 reads: 

^ = (1,1,-1,0), g 2 = (2,2,0,-1). 

Corresponding Eq. (18) is: 

A*h + A*o — /^OH, 2/iH + 2/io = A t H 2 o 2 - 

General construction of the quasi-equilibrium manifold: In the space of concentration, one 
defines a subspace L which satisfies the balance constraints: 

{KL) = 0. 

The orthogonal complement of L in the space with coordinates \i = VG defines then the 
quasi-equilibrium manifold Ql- For the actual computations, one requires the inversion 
from n to c. Duality structure fi <-> c is well studied by many authors [62,19]. 
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Quasi- equilibrium projector. It is not sufficient to just derive the manifold, it is also re- 
quired to define a projector which would transform the vector field defined on the space 
of concentrations to a vector field on the manifold. Quasi-equilibrium manifold consists 
of points which minimize G on the affine spaces of the form c + L. These affine planes 
are hypothetic planes of fast motions (G is decreasing in the course of the fast motions). 
Therefore, the quasi-equilibrium projector maps the whole space of concentrations on Ql 
parallel to L. The vector field is also projected onto the tangent space of Q,l parallel to 
L. 

Thus, the quasi-equilibrium approximation implies the decomposition of motions into the 
fast - parallel to L, and the slow - along the quasi-equilibrium manifold. In order to 
construct the quasi-equilibrium approximation, knowledge of reaction rate constants of 
"fast" reactions is not required (stoichiometric vectors of all these fast reaction are in L, 
Tfast e L, thus, knowledge of L suffices), one only needs some confidence in that they all 
are sufficiently fast [74] . The quasi-equilibrium manifold itself is constructed based on the 
knowledge of L and of G. Dynamics on the quasi-equilibrium manifold is defined as the 
quasi-equilibrium projection of the "slow component" of kinetic equations (4). 

2.4 Model equations 

The assumption behind the quasi-equilibrium is the hypothesis of the decomposition of 
motions into fast and slow. The quasi-equilibrium approximation itself describes slow 
motions. However, sometimes it becomes necessary to restore to the whole system, and 
to take into account the fast motions as well. With this, it is desirable to keep intact one 
of the important advantages of the quasi-equilibrium approximation - its independence of 
the rate constants of fast reactions. For this purpose, the detailed fast kinetics is replaced 
by a model equation (single relaxation time approximation) . 

Quasi-equilibrium models (QEM) are constructed as follows: For each concentration vector 
c, consider the affine manifold, c + L. Its intersection with the quasi-equilibrium manifold 
Ql consists of one point. This point delivers the minimum to G on c + L. Let us denote 
this point as c* L (c). The equation of the quasi-equilibrium model reads: 

c=--[c~ c* L (c)} + £ 7s W s (c* L (c)), (19) 

slow 

where r > is the relaxation time of the fast subsystem. Rates of slow reactions are 
computed in the points c* L (c) (the second term in the right hand side of Eq. (19), whereas 
the rapid motion is taken into account by a simple relaxational term (the first term in the 
right hand side of Eq. (19). The most famous model kinetic equation is the BGK equation 
in the theory of the Boltzmann equation [8]. The general theory of the quasi-equilibrium 
models, including proofs of their thermodynamic consistency, was constructed in [27,29]. 

Single relaxation time gradient models (SRTGM) were considered in [2,3,4] in the context 
of the lattice Boltzmann method for hydrodynamics. These models are aimed at improving 
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the obvious drawback of quasi-equilibrium models (19): In order to construct the QEM, 
one needs to compute the function, 

cUc) = arg min G(x). (20) 
LV ' ^xec+L, x>o v ' v ' 



This is a convex programming problem. It does not always has a closed-form solution. 

Let g 1: . . . ,g k is the orthonormal basis of L. We denote as D(c) the k x k matrix with 
the elements (g i: H c gj), where H c is the matrix of second derivatives of G (7). Let C(c) 
be the inverse of D(c). The single relaxation time gradient model has the form: 

c = -WoMcUgj, VG) + £ lsW s {c). (21) 

i,j slow 



The first term drives the system to the minimum of G on c + L, it does not require 
solving the problem (20), and its spectrum in the quasi-equilibrium is the same as in the 
quasi-equilibrium model (19). Note that the slow component is evaluated in the "current" 
state c. 

The models (19) and (21) lift the quasi-equilibrium approximation to a kinetic equation 
by approximating the fast dynamics with a single "reaction rate constant" - relaxation 
time r. 



2.5 Quasi-steady state approximation 



The quasi-steady state approximation (QSS) is a tool used in a huge amount of works. 
Let us split the list of species in two groups: The basic and the intermediate (radicals 
etc). Concentration vectors are denoted accordingly, c s (slow, basic species), and c f (fast, 
intermediate species). The concentration vector c is the direct sum, c = c s © c f . The fast 
subsystem is Eq. (4) for the component c f at fixed values of c s . If it happens that this way 
defined fast subsystem relaxes to a stationary state, c f — > c (cf), then the assumption 
that c f = Cq SS (c) is precisely the QSS assumption. The slow subsystem is the part of 
the system (4) for c s , in the right hand side of which the component c f is replaced with 
Cq SS (c). Thus, J = J s (B J{, where 



c f = J f (c s © c f ), c s = const; c f -> q L(c s ); (22) 



-'qss 

• S T ( « s ^ f 



c s = J s (c s ©c ss (c s )). (23) 



Bifurcations in the system (22) under variation of c s as a parameter are confronted to 
kinetic critical phenomena. Studies of more complicated dynamic phenomena in the fast 
subsystem (22) require various techniques of averaging, stability analysis of the averaged 
quantities etc. 
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Various versions of the QSS method are well possible, and are actually used widely, for 
example, the hierarchical QSS method. There, one defines not a single fast subsystem 
but a hierarchy of them, c fl , . . . , c ik . Each subsystem c il is regarded as a slow system 
for all the foregoing subsystems, and it is regarded as a fast subsystem for the following 
members of the hierarchy. Instead of one system of equations (22), a hierarchy of systems 
of lower- dimensional equations is considered, each of these subsystem is easier to study 
analytically. 

Theory of singularly perturbed systems of ordinary differential equations is used to provide 
a mathematical background and further development of the QSS approximation [10,68]. 
In spite of a broad literature on this subject, it remains, in general, unclear, what is the 
smallness parameter that separates the intermediate (fast) species from the basic (slow). 
Reaction rate constants cannot be such a parameter (unlike in the case of the quasi- 
equilibrium) . Indeed, intermediate species participate in the same reactions, as the basic 
species (for example, H 2 ^ 2H, H + 2 ^ OH + 0). It is therefore incorrect to state that 
c f evolve faster than c s . In the sense of reaction rate constants, c f is not faster. 

For catalytic reactions, it is not difficult to figure out what is the smallness parameter that 
separates the intermediate species from the basic, and which allows to upgrade the QSS 
assumption to a singular perturbation theory rigorously [77]. This smallness parameter 
is the ratio of balances: Intermediate species include the catalyst, and their total amount 
is simply significantly less than the amount of all the q's. After renormalizing to the 
variables of one order of magnitude, the small parameter appears explicitly. 

For usual radicals, the origin of the smallness parameter is quite similar. There are much 
less radicals than the basic species (otherwise, the QSS assumption is inapplicable). In 
the case of radicals, however, the smallness parameter cannot be extracted directly from 
balances Bi (2). Instead, one can come up with a thermodynamic estimate: Function 
G decreases in the course of reactions, whereupon we obtain the limiting estimate of 
concentrations of any specie: 

Ci < max &, (24) 

G(C)<G(C(0)) 

where c(0) is the initial composition. If the concentration cr of the radical R is small 
both initially and in the equilibrium, then it should remain small also along the path to 
the equilibrium. For example, in the case of ideal G (9) under relevant conditions, for any 
t > 0, the following inequality is valid: 

c R [ln(c R (t)/4 q ) - 1] < G(c(0)). (25) 

Inequality (25) provides the simplest (but rather coarse) thermodynamic estimate of c^(t) 
in terms of G(c(0)) and uniformly for t > 0. Complete theory of thermodynamic 
estimates of dynamics has been developed in [23]. One can also do computations without 
a priori estimations, if one accepts the QSS assumption until the values c f stay sufficiently 
small. 
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Let us assume that an a priori estimate has been found, < q max , for each q. 

These estimate may depend on the initial conditions, thermodynamic data etc. With 
these estimates, we are able to renormalize the variables in the kinetic equations (4) in 
such a way that renormalized variables take their values from the unit segment [0,1]: 
q = Ci/ci max . Then the system (4) can be written as follows: 



The system of dimensionless parameters, = q max / maxj q max defines a hierarchy of 
relaxation times, and with its help one can establish various realizations of the QSS ap- 
proximation. The simplest version is the standard QSS assumption: Parameters are 
separated in two groups, the smaller ones, and of the order 1. Accordingly, the concen- 
tration vector is split into c s © c f . Various hierarchical QSS are possible, with this, the 
problem becomes more tractable analytically. 

Corrections to the QSS approximation can be addressed in various ways (see, e. g., [72,70]). 
There exist a variety of ways to introduce the smallness parameter into kinetic equations, 
and one can find applications to each of the realizations. However, the two particular 
realizations remain basic for chemical kinetics: (i) Fast reactions (under a given ther- 
modynamic data); (ii) Small concentrations. In the first case, one is led to the quasi- 
equilibrium approximation, in the second case - to the classical QSS assumption. Both 
of these approximations allow for hierarchical realizations, those which include not just 
two but many relaxation time scales. Such a multi-scale approach essentially simplifies 
analytical studies of the problem. 

The method of invariant manifold which we present below in the section 6 allows to use 
both the QE and the QSS as initial approximations in the iterational process of seeking 
slow invariant manifolds. It is also possible to use a different initial ansatz chosen by a 
physical intuition, like, for example, the Tamm-Mott- Smith approximation in the theory 
of strong shock waves [24]. 

2.6 Methods based on spectral decomposition of Jacobian fields 

The idea to use the spectral decomposition of Jacobian fields in the problem of separating 
the motions into fast and slow originates from methods of analysis of stiff systems [22], and 
from methods of sensitivity analysis in control theory [64] . There are two basic statements 
of the problem for these methods: (i) The problem of the slow manifold, and (ii) The 
problem of a complete decomposition (complete integrability) of kinetic equations. The 
first of these problems consists in constructing the slow manifold fi, and a decomposition 
of motions into the fast one - towards S7, and the slow one - along Q [61]. The second of 
these problems consists in a transformation of kinetic equations (4) to a diagonal form, 
d = fi(d) (so-called full nonlinear lumping or modes decoupling, [56,59,71]). Clearly, if 
one finds a sufficiently explicit solution to the second problem, then the system (4) is 
completely integrable, and nothing more is needed, the result has to be simply used. The 




do. 



1 



(26) 



15 



question is only to what extend such a solution can be possible, and how difficult it would 
be as compared to the first problem to find it. 

One of the currently most popular methods is the construction of the so-called intrinsic 
low- dimensional manifold (ILDM, [61]). This method is based on the following geometric 
picture: For each point c, one defines the Jacobian matrix of Eq. (4), F c = dJ(c)/dc. 
One assumes that, in the domain of interest, the eigenvalues of Fc are separated into two 
groups, and A^, and that the following inequalities are valid: 

Re A ■ > a > 6 > ReAj, a > b, b < 0. 

Let us denote as L s c and L c the invariant subspaces corresponding to A s and A f , re- 
spectively, and let Z s c and Z c be the corresponding spectral projectors, Z S C L S C = L s c , 
Z* C ]J C = L f c , Z s c L f c = Z f c L s c = {0}, Z s c + Z f c = 1. Operator Z s c projects onto the 
subspace of "slow modes" L s c , and it annihilates the "fast modes" L { c . Operator Z c does 
the opposite, it projects onto fast modes, and it annihilates the slow modes. The basic 
equation of the ILDM reads: 

Z f c J(c) = 0. (27) 

In this equation, the unknown is the concentration vector c. The set of solutions to Eq. 
(27) is the ILDM manifold fi ildm . 

For linear systems, F c , Z s c , and Z C) do not depend on c, and flwdm = c cq + L s . On the 
other hand, obviously, c cq e findm- Therefore, procedures of solving of Eq. (27) can be 
initiated by choosing the linear approximation, f2||^ m = c cq + L s C c q , in the neighborhood 
of the equilibrium c cq , and then continued parametrically into the nonlinear domain. 
Computational technologies of a continuation of solutions with respect to parameters are 
well developed (see, for example, [55,65]). The problem of the relevant parameterization 
is solved locally: In the neighborhood of a given point c° one can choose Z s c (c — c°) 
for a characterization of the vector c. In this case, the space of parameters is L s c . There 
exist other, physically motivated ways to parameterize manifolds ([24]; see also section 
4.1 below). 

There are two drawbacks of the ILDM method which call for its refinement: (i) "Intrin- 
sic" does not imply "invariant" . Eq. (27) is not invariant of the dynamics (4). If one 
differentiates Eq. (27) in time due to Eq. (4), one obtains a new equation which is the 
implication of Eq. (27) only for linear systems. In a general case, the motion c(t) takes off 
the fiiidm- Invariance of a manifold fl means that J(c) touches Q in every point c e f2. 
It remains unclear how the ILDM (27) corresponds with this condition. Thus, from the 
dynamical perspective, the status of the ILDM remains not well defined, or "ILDM is 
ILDM", defined self-consistently by Eq. (27), and that is all what can be said about it. 
(ii) From the geometrical standpoint, spectral decomposition of Jacobian fields is not the 
most attractive way to compute manifolds. If we are interested in the behavior of trajec- 
tories, how they converge or diverge, then one should consider the symmetrized part of 
Fc, rather than Fc itself. 



16 



Symmetric part, F^ 111 = (1/2) (F* c + Fq), defines the dynamics of the distance between 
two solutions, c and c', in a given local Euclidean metrics. Skew-symmetric part de- 
fines rotations. If we want to study manifolds based on the argument about conver- 
gence/divergence of trajectories, then we should use in Eq. (27) the spectral projector 
Zc ym for the operator F^™ 1 . This, by the way, is also a significant simplification from the 
standpoint of computations. It remains to choose the metrics. This choice is unambiguous 
from the thermodynamic perspective. In fact, there is only one choice which fits into the 
physical meaning of the problem, this is the metrics associated with the thermodynamic 
(or entropic) scalar product, 

(x,y) = (x,H c y), (28) 



where He is the matrix of second-order derivatives of G (7). In the equilibrium, operator 
Fc cc < is selfadjoint with respect to this scalar product (Onsager's reciprocity relations). 
Therefore, the behavior of the ILDM in the vicinity of the equilibrium does not alter 
under the replacement, F c ^ = F s ^. In terms of usual matrix representation, we have: 

F s r = l(Fc + H c 'F T c H c ), (29) 



where F T C is the ordinary transposition. 

The ILDM constructed with the help of the symmetrized Jacobian field will be termed 
the symmetric entropic intrinsic low- dimensional manifold (SEILDM). Selfadjointness of 
F s ^ m (29) with respect to the thermodynamic scalar product (28) simplifies considerably 
computations of spectral decomposition. Moreover, it becomes necessary to do spectral 
decomposition in only one point - in the equilibrium. Perturbation theory for selfadjoint 
operators is a very well developed subject [53], which makes it possible to easily extend 
the spectral decomposition with respect to parameters. A more detailed discussion of the 
selfadjoint linearization will be given below in section 5.2. 

Thus, when the geometric picture behind the decomposition of motions is specified, the 
physical significance of the ILDM becomes more transparent, and it leads to its modifi- 
cation into the SEILDM. This also gains simplicity in the implementation by switching 
from non-self adjoint spectral problems to selfadjoint. The quantitative estimate of this 
simplification is readily available: Let d be the dimension of the phase space, and k the 
dimension of the ILDM (k = dimL^,). The space of all the projectors Z with the k- 
dimensional image has the dimension D = 2k(d — k). The space of all the selfadjoint 
projectors with the k- dimensional image has the dimension D sym = k(d — k). For d = 20 
and k = 3, we have D = 102 and D sym = 51. When the spectral decomposition by means 
of parametric extension is addressed, one considers equations of the form: 

^ = W£, Z iM ,F cw ,VF cw ), (30, 



where r is the parameter, and V-Fc = VVJ(c) is the differential of the Jacobian field. 
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For the selfadjoint case, where we use = F s ^ m instead of Fc, this system of equations has 
twice less independent variables, and also the right hand is of a simpler structure. 

It is more difficult to improve on the first of the remarks (ILDM is not invariant). The 
following naive approach may seem possible: 

(i) Take f2 ildm = c cq + L s C eq in a neighborhood U of the equilibrium c cq . [This is also a 
useful initial approximation for solving Eq. (27)]. 

(ii) Instead of computing the solution to Eq. (27), integrate the kinetic equations (4) 
backwards in the time. It is sufficient to take initial conditions c(0) from a dense set on 
the boundary, dU n (c cq + L^eq), and to compute solutions c(t), t < 0. 

(iii) Consider the obtained set of trajectories as an approximation of the slow invariant 
manifold. 

This approach will guarantee invariance, by construction, but it is prone to pitfalls in what 
concerns the slowness. Indeed, the integration backwards in the time will see exponentially 
divergent trajectories, if they were exponentially converging in the normal time progress. 
This way one finds some invariant manifold which touches c cq + L s C c q in the equilibrium. 
Unfortunately, there are infinitely many such manifolds, and they fill out almost all the 
space of concentrations. However, we must select the slow component of motions. Such 
a regularization is possible. Indeed, let us replace in Eq. (4) the vector field J(c) by the 
vector field Zc ym J(c), and obtain a regularized kinetic equation, 

c = Z s * ym J(c). (31) 

Let us replace integration backwards in time of the kinetic equation (4) in the naive 
approach described above by integration backwards in time of the regularized kinetic 
equation (31). With this, we obtain a rather convincing version of the ILDM (SEILDM). 
Using Eq. (30), one also can write down an equation for the projector Z s ^ ym , putting r — t. 
Replacement of Eq. (4) by Eq. (31) also makes the integration backwards in time in the 
naive approach more stable. However, regularization will again conflict with invariance. 
The "naive refinement" after the regularization (31) produces just a slightly different ver- 
sion of the ILDM (or SEILDM) but it does not construct the slow invariant manifold. So, 
where is the way out? We believe that the ILDM and its version SEILDM are, in general, 
good initial approximations of the slow manifold. However, if one is indeed interested in 
finding the invariant manifold, one has to write out the true condition of invariance and 
solve it. As for the initial approximation for the method of invariant manifold one can use 
any ansatz, in particular, the SEILDM. 

The problem of a complete decomposition of kinetic equations can be solved indeed in some 
cases. The first such solution was the spectral decomposition for linear systems [75]. De- 
composition is sometimes possible also for nonlinear systems ([59]; [71]). The most famous 
example of a complete decomposition of infinite-dimensional kinetic equation is the com- 
plete integrability of the space-independent Boltzmann equation for Maxwell's molecules 
found in [9]. However, in a general case, there exist no analytical, not even a twice differen- 
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tiable transformation which would decouple modes. The well known Grobman-Hartman 
theorem [43,44] states only the existence of a continuous transform which decomposes 
modes in a neighborhood of the equilibrium. For example, the analytic planar system, 
dx/dt = —x, dy/dt = —2y + x 2 , is not C 2 linearizable. These problems remain of interest 
[15]. Therefore, in particular, it becomes quite ineffective to construct such a transforma- 
tion in a form of a series. It is more effective to solve a simpler problem of extraction of 
a slow invariant manifold [7] . 

Sensitivity analysis [64,63,56] makes it possible to select essential variables and reactions, 
and to decompose motions into fast and slow. In a sense, the ILDM method is a devel- 
opment of the sensitivity analysis. In particular, the computational singular perturbation 
(CSP) method of [56] includes ILDM (or any other reasonable initial choice of the mani- 
fold) into a procedure of consequent refinements. Recently, a further step in this direction 
was done in [78] . In this work, the authors use a nonlocal in time criterion of closeness 
of solutions of the full and of the reduced systems of chemical kinetics. They require not 
just a closeness of derivatives but a true closeness of the dynamics. 

Let us be interested in the dynamics of the concentrations of just a few species, A 1; . . . , A p , 
whereas the rest of the species, A p+1 , . . . , A n are used for building the kinetic equation, 
and for understanding the process. Let c goa i be the concentration vector with components 
Ci, . . . , c p , Cg 0a i(t) be the corresponding components of the solution to Eq. (4), and 
be the solution to the simplified model with corresponding initial conditions. [78] suggest 
to minimize the difference between c goa i(t) and on the segment t G [0, T\. ||c goa i(t) — 
c goaill ~~ * mm - I n the course of the optimization under certain restrictions one selects the 
optimal (or appropriate) reduced model. The sequential quadratic programming method 
and heuristic rules of sorting the reactions, substances etc were used. In the result, for 
some stiff systems studied, one avoids typical pitfalls of the local sensitivity analysis. In 
simpler situations this method should give similar results as the local methods. 



2. 7 Thermodynamic criteria for selection of important reactions 



One of the problems addressed by the sensitivity analysis is the selection of the important 
and discarding the unimportant reactions. [12] suggested a simple principle to compare 
importance of different reactions according to their contribution to the entropy production 
(or, which is the same, according to their contribution to G). Based on this principle, [17] 
described domains of parameters in which the reaction of hydrogen oxidation, H2 + O2 + M, 
proceeds due to different mechanisms. For each elementary reaction, he has derived the 
domain inside which the contribution of this reaction is essential (nonnegligible). Due to 
its simplicity, this entropy production principle is especially well suited for analysis of 
complex problems. In particular, recently, a version of the entropy production principle 
was used in the problem of selection of boundary conditions for Grad's moment equations 
[69,42]. For ideal systems (9), the contribution of the sth reaction to G has a particularly 
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simple form: 



G s = -W s In 



( 



w- 



) 



s=l 



r 



(32) 



For nonideal systems, the corresponding expressions (13) are also not too complicated. 



3 Outline of the method of invariant manifold 

In many cases, dynamics of the <i-dimensional system (4) leads to a manifold of a lower 
dimension. Intuitively, a typical phase trajectory behaves as follows: Given the initial 
state c(0) at t — 0, and after some period of time, the trajectory comes close to some 
low-dimensional manifold f2, and after that proceeds towards the equilibrium essentially 
along this manifold. The goal is to construct this manifold. 

The starting point of our approach is based on a formulation of the two main requirements: 

(i) . Dynamic invariance: The manifold O should be (positively) invariant under the dy- 
namics of the originating system (4): If c(0) G fi, then c(t) G Q for each t > 0. 

(ii) . Thermodynamic consistency of the reduced dynamics: Let some (not obligatory in- 
variant) manifold f2 is considered as a manifold of reduced description. We should define 
a set of linear operators, Pc, labeled by the states c G f2, which project the vectors J(c), 
c G Q onto the tangent bundle of the manifold f2, thereby generating the induced vector 
field, PcJ{c), c G CI. This induced vector field on the tangent bundle of the manifold 
O is identified with the reduced dynamics along the manifold fl. The thermodynamicity 
requirement for this induced vector field reads 



In order to meet these requirements, the method of invariant manifold suggests two com- 
plementary procedures: 

(i) . To treat the condition of dynamic invariance as an equation, and to solve it iteratively 
by a Newton method. This procedure is geometric in its nature, and it does not use the 
time dependence and small parameters. 

(ii) . Given an approximate manifold of reduced description, to construct the projector 
satisfying the condition (33) in a way which does not depend on the vector field J. 

We shall now outline both these procedures starting with the second. The solution consists, 
in the first place, in formulating the thermodynamic condition which should be met by 
the projectors P C '. For each c G £1, let us consider the linear functional 




(33) 




(34) 
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Then the thermodynamic condition for the projectors reads: 
kerP c C kerM£, for each c e ft. 



(35) 



Here kerPc is the null space of the projector, and kerM£ is the hyperplane orthogonal 
to the vector M* c . It has been shown [24,28] that the condition (35) is the necessary 
and sufficient condition to establish the thermodynamic induce vector field on the given 
manifold fl for all possible dissipative vector fields J simultaneously. 

Let us now turn to the requirement of invariance. By a definition, the manifold fl is 
invariant with respect to the vector field J if and only if the following equality is true: 

[1 - P] J(c) = 0, for each c e ft. (36) 

In this expression P is an arbitrary projector on the tangent bundle of the manifold ft. It 
has been suggested to consider the condition (36) as an equation to be solved iteratively 
starting with some appropriate initial manifold. 

Iterations for the invariance equation (36) are considered in the section 5. The next section 
presents construction of the thermodynamic projector using a specific parameterization 
of manifolds. 



4 Thermodynamic projector 

4-1 Thermodynamic parameterization 

In this section, ft denotes a generic p-dimensional manifold. First, it should be men- 
tioned that any parameterization of fl generates a certain projector, and thereby a 
certain reduced dynamics. Indeed, let us consider a set of m independent functionals 
M(c) = {Mi(c), . . . , M p (c)}, and let us assume that they form a coordinate system on 
f2 in such a way that fl = c(M), where c(M) is a vector function of the parameters 
Mi, . . . , M p . Then the projector associated with this parameterization reads: 

v Q C (M) 

Pc { m)X = E^^(VM, | c(M) ,aO, (37) 

i=l % 

where N^ 1 is the inverse to the p x p matrix: 

N(M) = IKVM^dc/dM^l (38) 

This somewhat involved notation is intended to stress that the projector (37) is dictated by 
the choice of the parameterization. Subsequently, the induced vector field of the reduced 
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dynamics is found by applying projectors (37) on the vectors J(c(M)), thereby inducing 
the reduced dynamics in terms of the parameters M as follows: 

M t = £ iV i - 1 (M)(VM, \ C(M) , J(c(M))), (39) 
j'=i 



Depending on the choice of the parameterization, dynamic equations (39) are (or are not) 
consistent with the thermodynamic requirement (33). The thermodynamic parameteriza- 
tion makes use of the condition (35) in order to establish the thermodynamic projector. 
Specializing to the case (37), let us consider the linear functionals, 

DM t \ ciM) (x) = (VM l \ c(M) ,x). (40) 



Then the condition (35) takes the form: 

HkerDM, | c(M) CkerM* (M) , (41) 



i=i 



that is, the intersection of null spaces of the functionals (40) should belong to the null 
space of the differential of the Lyapunov function G, in each point of the manifold f2. 

In practice, in order to construct the thermodynamic parameterization, we take the fol- 
lowing set of functionals in each point c of the manifold f2: 



M 1 (x) = M* c (x), ceCl (42) 
M i (x) = (m i ,x), i = 2,...,p (43) 

It is required that vectors VG(c), ra 2 , . . . , m p are linearly independent in each state 
c G fi. Inclusion of the functionals (34) as a part of the system (42) and (43) implies the 
thermodynamic condition (41). Also, any linear combination of the parameter set (42), 
(43) will meet the thermodynamicity requirement. 

It is important to notice here that the thermodynamic condition is satisfied whatsoever 
the functionals M 2 , . . . , M p are. This is very convenient for it gives an opportunity to take 
into account the conserved quantities correctly. The manifolds we are going to deal with 
should be consistent with the conservation laws (2). While the explicit characterization 
of the phase space V is a problem on its own, in practice, it is customary to work in the 
n-dimensional space while keeping the constraints (2) explicitly on each step of the con- 
struction. For this technical reason, it is convenient to consider manifolds of the dimension 
p > I, where I is the number of conservation laws, in the n-dimensional space rather than 
in the phase space V. The thermodynamic parameterization is then concordant also with 
the conservation laws if / of the linear functionals (43) are identified with the conservation 
laws. In the sequel, only projectors consistent with conservation laws are considered. 
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Very frequently, the manifold fl is represented as a p-parametric family c(ai, . . . , a p ), 
where are coordinates on the manifold. The thermodynamic re-parameterization sug- 
gests a representation of the coordinates in terms of M c , M 2 , . . . , M p (42), (43). While 
the explicit construction of these functions may be a formidable task, we notice that the 
construction of the thermodynamic projector of the form (37) and of the dynamic equa- 
tions (39) is relatively easy because only the derivatives dc/dMi enter these expressions. 
This point was discussed in a detail in [24,28]. 

4-2 Decomposition of motions: Thermodynamics 

Finally, let us discuss how the thermodynamic projector is related to the decomposition 
of motions. Assuming that the decomposition of motions near the manifold fl is true 
indeed, let us consider states which were initially close enough to the manifold ft. Even 
without knowing the details about the evolution of the states towards O, we know that 
the Lyapunov function G was decreasing in the course of this evolution. Let us consider 
a set of states Uc which contains all those vectors c! that have arrived (in other words, 
have been projected) into the point c e O. Then we observe that the state c furnishes 
the minimum of the function G on the set Uc- If a state c' e Uc-, and if it deviates small 
enough from the state c so that the linear approximation is valid, then c' belongs to the 
affine hyperplane 

T c = c + ker M* , c E CI. (44) 

This hyperplane actually participates in the condition (35). The consideration was entitled 
'thermodynamic' [24] because it describes the states c G f2 as points of minimum of the 
function G over the corresponding hyperplanes (44). 



5 Corrections 

5.1 Preliminary discussion 

The thermodynamic projector is needed to induce the dynamics on a given manifold 
in such a way that the dissipation inequality (33) holds. Coming back to the issue of 
constructing corrections, we should stress that the projector participating in the invariance 
condition (36) is arbitrary. It is convenient to make use of this point: When Eq. (36) 
is solved iteratively, the projector may be kept non-thermodynamic unless the induced 
dynamics is explicitly needed. 

Let us assume that we have chosen the initial manifold, Oo, together with the associated 
projector P , as the first approximation to the desired manifold of reduced description. 
Though the choice of the initial approximation f2 depends on the specific problem, it is 



23 



often reasonable to consider quasi-equilibrium or quasi steady-state approximations. In 
most cases, the manifold O is not an invariant manifold. This means that Qq does not 
satisfy the invariance condition (36): 



A = [1 — Pq] J(cq) 7^ 0, for some c G fto- 



(45) 



Therefore, we seek a correction C\ = c + 5c. Substituting P = P and c = c + Sc into 
the invariance equation (36), and after the linearization in Sc, we derive the following 
linear equation: 



where L Co is the matrix of first derivatives of the vector function J, computed in the 
state c G f2o- The system of linear algebraic equations (46) should be supplied with the 
additional condition. 



In order to illustrate the nature of the Eq. (46), let us consider the case of linear manifolds 
for linear systems. Let a linear evolution equation is given in the finite-dimensional real 
space: c = Lc, where L is negatively definite symmetric matrix with a simple spectrum. 
Let us further assume the quadratic Lyapunov function, G(c) = (c, c). The manifolds 
we consider are lines, 1(a) = ae, where e is the unit vector, and a is a scalar. The 
invariance equation for such manifolds reads: e(e, Le) — Le = 0, and is simply a form of 
the eigenvalue problem for the operator L. Solutions to the latter equation are eigenvectors 
ej, corresponding to eigenvalues Aj. 

Assume that we have chosen a line, Iq = ae , defined by the unit vector e , and that eo is 
not an eigenvector of L. We seek another line, l\ = aei, where ei is another unit vector, 
e i — 2/i/l|2/ill) Vi — e o + $y- The additional condition (47) now reads: (5y,e ) = 0. Then 
the Eq. (46) becomes [1 — e (e , -)]-M e o + Sy] = 0. Subject to the additional condition, 
the unique solution is as follows: eo + 5y = (eo, L~ l eo)~ 1 L~ 1 eo. Rewriting the latter 
expression in the eigen-basis of L, we have: e + 5y oc J2i \ rle i( e ii e o)- The leading term 
in this sum corresponds to the eigenvalue with the minimal absolute value. The example 
indicates that the method of linearization (46) seeks the direction of the slowest relaxation. 
For this reason, the method (46) can be recognized as the basis of an iterative method 
for constructing the manifolds of slow motions. 

For the nonlinear systems, the matrix L Co in the Eq. (46) depends nontrivially on Cq. In 
this case the system (46) requires a further specification which will be done now. 



[1 - P ] [J (cq) + LcoSc] = 0, 



(46) 



P Q 5c = 0. 



(47) 
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5.2 Symmetric linearization 



The invariance condition (36) supports a lot of invariant manifolds, and not all of them 
are relevant to the reduced description (for example, any individual trajectory is itself an 
invariant manifold). This should be carefully taken into account when deriving a relevant 
equation for the correction in the states of the initial manifold Oo which are located far 
from equilibrium. This point concerns the procedure of the linearization of the vector field 
J, appearing in the equation (46). We shall return to the explicit form of the Marcelin-De 
Donder kinetic function (11). Let c is an arbitrary fixed element of the phase space. The 
linearization of the vector function J (12) about c may be written J(c+Sc) ps J(c)+LqSc 
where the linear operator Lc acts as follows: 

L c x = j^ ls [W:{c){cx s ,H c x) - W-(c)(f3 s ,H c x)]. (48) 

s=l 



Here H c is the matrix of second derivatives of the function G in the state c [see Eq. (7)]. 
The matrix L c in the Eq. (48) can be decomposed as follows: 

L C = L' C + L' C . (49) 



Matrices L' c and L" c act as follows: 



L 'c* = ~\ + W7(c)] 7 .(7„ **cx), (50) 

Z s=l 

L'c* = \ E[^ s + (c) - W-(c)] ls ( as + (3 S , H c x). (51) 

Z 8=1 

Some features of this decomposition are best seen when we use the thermodynamic scalar 
product (28): The following properties of the matrix L' c are verified immediately: 

(i) The matrix L' c is symmetric in the scalar product (28): 

(x,L' c y) = (y,L' c x). (52) 



(ii) The matrix L' c is nonpositive definite in the scalar product (28): 

(x,L' c x)<0. (53) 

(iii) The null space of the matrix L' c is the linear envelope of the vectors H^bi repre- 
senting the complete system of conservation laws: 

keri' c = LmiH^bi, i = 1, . . . , 1} (54) 
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(iv) If c = c°^, then W+(c cq ) = W-(c eq ), and 




(55) 



Thus, the decomposition Eq. (49) splits the matrix L c in two parts: one part, Eq. (50) 
is symmetric and nonpositive definite, while the other part, Eq. (51), vanishes in the 
equilibrium. The decomposition Eq. (49) explicitly takes into account the Marcelin-Dc 
Donder form of the kinetic function. For other dissipative systems, the decomposition 
(49) is possible as soon as the relevant kinetic operator is written in a gain-loss form [for 
instance, this is straightforward for the Boltzmann collision operator]. 

In the sequel, we shall make use of the properties of the operator L' c (50) for constructing 
the dynamic correction by extending the picture of the decomposition of motions. 



5.3 Decomposition of motions: Kinetics 

The assumption about the existence of the decomposition of motions near the manifold 
of reduced description ft has led to the thermodynamic specifications of the states c G 
O. This was accomplished in the section 4.2, where the thermodynamic projector was 
backed by an appropriate variational formulation, and this helped us to establish the 
induced dynamics consistent with the dissipation property. Another important feature of 
the decomposition of motions is that the states c G O can be specified kinetically. Indeed, 
let us do it again as if the decomposition of motions were valid in the neighborhood of the 
manifold ft, and let us 'freeze' the slow dynamics along the ft, focusing on the fast process 
of relaxation towards a state c G f2. From the thermodynamic perspective, fast motions 
take place on the afline hyperplane c + 5c G r Co , where r C() is given by Eq. (44). From 
the kinetic perspective, fast motions on this hyperplane should be treated as a relaxation 
equation, equipped with the quadratic Lyapunov function 5G = (5c, 5c), Furthermore, 
we require that the linear operator of this evolution equation should respect Onsager's 
symmetry requirements (selfadjointness with respect to the entropic scalar product). This 
latter crucial requirement describes fast motions under the frozen slow evolution in the 
similar way, as all the motions near the equilibrium. 

Let us consider now the manifold Qq which is not the invariant manifold of the reduced 
description but, by our assumption, is located close to it. Consider a state c G f2 , and 
the states c + 5c close to it. Further, let us consider an equation 



Due to the properties of the operator L' C( (50), this equation can be regarded as a model 
of the assumed true relaxation equation near the true manifold of the reduced description. 
For this reason, we shall use the symmetric operator L' c (50) instead of the linear operator 
Lc when constructing the corrections. 



5 c = L' r 5c. 



(56) 
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5.4 Symmetric iteration 



Let the manifold f2 and the corresponding projector P are the initial approximation to 
the invariant manifold of the reduced description. The dynamic correction C\ = c + 5c 
is found upon solving the following system of linear algebraic equations: 

= 0, P 5c = 0. (57) 

Here L' Co is the matrix (50) taken in the states on the manifold fio- An important technical 
point here is that the linear system (57) always has the unique solution for any choice of 
the manifold f2. This point is crucial since it guarantees the opportunity of carrying out 
the correction process for arbitrary number of steps. 



6 The method of invariant manifold 

We shall now combine together the two procedures discussed above. The resulting method 
of invariant manifold intends to seek iteratively the reduced description, starting with an 
initial approximation. 

(i). Initialization. In order to start the procedure, it is required to choose the initial 
manifold fi , and to derive corresponding thermodynamic projector P . In the majority 
of cases, initial manifolds are available in two different ways. The first case are the quasi- 
equilibrium manifolds described in the section 2.3. The macroscopic parameters are Mj = 
q = (mj,c), where rrii is the unit vector corresponding to the specie A{. The quasi- 
equilibrium manifold, c (M 1? . . . , M k , B ± , . . . , B t ), compatible with the conservation laws, 
is the solution to the variational problem: 

G — > min , (m^, c) = q, i — 1, . . . , k, (58) 
(bj,c) = Bj, j = 1, . . . ,1. 

In the case of quasi-equilibrium approximation, the corresponding thermodynamic pro- 
jector can be written most straightforwardly in terms of the variables Mf. 

Po* = t^(rn u x) + j:^-(b l , X ). (59) 

i=l % i=l % 

For quasi-equilibrium manifolds, a reparameterization with the set (42), (43) is not nec- 
essary ([24]; [28]). 

The second source of initial approximations are quasi-stationary manifolds (section 2.5). 
Unlike the quasi-equilibrium case, the quasi- stationary manifolds must be reparameterized 
in order to construct the thermodynamic projector. 



[1-Po] J(c ) + L'5c 
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(ii). Corrections. Iterations are organized in accord with the rule: If c m is the mth ap- 
proximation to the invariant manifold, then the correction c m+1 = c m + 5c is found from 
the linear algebraic equations, 

[l-P m ](J(c m ) + L' c Jc)=0, (60) 
P m 5c = 0. (61) 

Here L' Cm is the symmetric matrix (50) evaluated at the mth approximation. The projector 
P m is not obligatory thermodynamic at that step, and it is taken as follows: 

P m x = Y,^(rn t ,x)+Y / ^(b h x). (62) 

i=l 1 i=l % 



(iii). Dynamics. Dynamics on the mth manifold is obtained with the thermodynamic 
re-parameterization. 

In the next section we shall illustrate how this all works. 



7 Illustration: Two-step catalytic reaction 

Here we consider a two-step four- component reaction with one catalyst A 2 : 

A 1 + A 2 ^ A 3 ^ A 2 + A 4 . (63) 

We assume the Lyapunov function of the form (9), G — Y^t=i Q[l n (Q/Ci q ) — 1]. The kinetic 
equation for the four-component vector of concentrations, c = (ci, c 2 , c 3 , c 4 ), has the form 

c = ll W 1 + j 2 W 2 . (64) 

Here j 12 are stoichiometric vectors, 

7l = (-1,-1,1,0), 72 = (0,1,-1,1), (65) 

while functions Wi t2 are reaction rates: 

Wi = kfcic 2 - k^c 3 , W 2 = k^c-i - k^c 2 c A . (66) 

Here kf 2 are reaction rate constants. The system under consideration has two conservation 
laws, 

c 1 + c 3 + c 4 = B ± , c 2 + c 3 = B 2 , (67) 
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or (61,2, c) = £?i )2 , where 61 = (1, 0, 1, 1) and b 2 = (0, 1,1,0). The nonlinear system (64) 
is effectively two-dimensional, and we consider a one-dimensional reduced description. 

We have chosen the concentration of the specie A 1 as the variable of reduced description: 
M = Ci, and c\ = (m, c), where m = (1,0,0,0). The initial manifold c (M) was taken 
as the quasi-equilibrium approximation, i.e. the vector function c is the solution to the 
problem: 

G — > min for (m, c) = c 1: (b l5 c) = B 1: (b 2 , c) = B 2 . (68) 



The solution to the problem (68) reads: 



c 0l = d, (69) 
c 02 = B 2 - 0(ci), 

C03 = 0(Cl), 

0)4 = ^1 - Ci - 0(ci), 

0(M) = A(ci) - y^cO-B^-d), 

, g 2 (g 1 -c? > )+^( C r'+^- Cl ) 

The thermodynamic projector associated with the manifold (69) reads: 

P x = ^(m, x) + f|(&i, *) + ||(6 2 , x). (70) 



Computing A = [1 — P ]J(c ) we find that the inequality (45) takes place, and thus the 
manifold c is not invariant. The first correction, c\ = c + 5c, is found from the linear 
algebraic system (60) 



(l-P )L' 6c=-[l-P ]J(co), (71) 

5d = 
5ci + 5c 3 + 5c 4 = 

Sc 3 + 5c 2 = 0, (72) 

where the symmetric 4x4 matrix L' Q has the form (we write instead of Cq in the subscript 
in order to simplify notations): 

W^(co) + Wr(co) 7ll W 2 + (c ) + W 2 (c ) l2l 
2 c ; 2 c / 



The explicit solution Ci(ci,Bi,B 2 ) to the linear system (71) is easily found, and we do 
not reproduce it here. The process was iterated. On the k + 1 iteration, the following 
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projector P k was used: 



P k x = 



dc k 



(ra, x) + 



dc k 



(bi,x) + 



dc k 



(b 2 ,x). 



(74) 



dci 



dB l 



dB 2 



Notice that projector P k (74) is the thermodynamic projector only if k = 0. As we have 
already mentioned it above, in the process of finding the corrections to the manifold, the 
non-thermodynamic projectors are allowed. The linear equation at the k + 1 iteration is 
thus obtained by replacing c , Pq, and L' Q with c k , P k , and L' k in all the entries of the 



Once the manifold c k was obtained on the kth iteration, we derived the corresponding 
dynamics by introducing the thermodynamic parameterization (and the corresponding 
thermodynamic projector) with the help of the function (42). The resulting dynamic 
equation for the variable C\ in the /cth approximation has the form: 



Analytic results were compared with the results of the numerical integration. The following 
set of parameters was used: 



Direct numerical integration of the system has demonstrated that the manifold C3 = Cg q 
in the plane (ci,c^) attracts all individual trajectories. Thus, the reduced description in 
this example should extract this manifold. 

Fig. 1 demonstrates the quasi-equilibrium manifold (69) and the first two corrections 
found analytically. It should be stressed that we spend no special effort on the construc- 
tion of the initial approximation, that is, of the quasi-equilibrium manifold, have not used 
any information about the Jacobian field (unlike, for example, the ILDM or CSP methods 
discussed above) etc. It is therefore not surprising that in this way chosen initial quasi- 
equilibrium approximation is in a rather poor agreement with the reduced description. 
However, it should be appreciated that the further corrections rapidly improve the situa- 
tion while no small parameter considerations were used. This confirms our expectation of 
the advantage of using the iteration methods in comparison to methods based on a small 
parameter expansions for model reduction problems. 



Eqs. (71) and (73). 




(75) 



Here [VG \ Ck } t = ln[c kl /^]. 




0.4, k^ = 1.0; 
0.1, cf = 0.4, 
1.0, B 2 = 0.2. 
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8 Relaxation methods 



Relaxation method is an alternative to the Newton iteration method described in section 
5. It is a one-dimensional Galerkin approximation for the linearized invariance equation 
(46,47). We shall solve the invariance equation (46,47) (or symmetric invariance equation 
(60,61) in projection on the defect of invariance (45) A = [1 — P c ]J(c). 

Let f2 be the initial approximation to the invariant manifold, and we seek the first 
correction, C\ = c + Ti(co)A(co), where function t(cq) has a dimension of the time, and 
is found from the condition that the linearized vector field attached to the points of the 
new manifold is orthogonal to the initial defect, 



Further steps r k (c) are found in the same way. It is clear from the latter equations that 
the step of the relaxation method is equivalent to the Galerkin approximation for solving 
the step of the Newton method. Actually, the relaxation method was first introduced in 
these terms in [52]. An advantage of equation (77) is the explicit form of the size of the 
steps r k (c). This method was successfully applied to the Fokker-Plank equation [52]. 



9 Method of invariant manifold without a priori parameterization 

Formally, the method of invariant manifold does not require a global parameterization of 
the manifolds. However, in most of the cases, one makes use of a priori defined "macro- 
scopic" variables M. This is motivated by the choice of quasi-equilibrium initial approxi- 
mations. 

Let a manifold Q, be defined in the phase space of the system, its tangent space in the 
point c be T c fi. How to define the projector of the whole concentrations space onto T c fi 
without using any a priori parameterization of ft? 

The basis of the answer to this question is the condition of thermodynamicity (35). Let 
us denote E as the concentration space, and consider the problem of the choice of the 
projector in the quadratic approximation to the thermodynamic potential G: 



(A(c ), (1 - P Co )[J(c ) + r 1 (co)( J DcJ)c A C() ]) Co = 0. 



(76) 



Explicitly, 




(77) 



G q = (g, H c Ac) + I (Ac H c Ac) = (g, Ac) + I (Ac, Ac), 



(78) 
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where He is the matrix of the second-order derivatives of G (7), g — iJ c 1 VG, Ac is the 
deviation of the concentration vector from the expansion point. 

Let a linear subspace T be given in the concentrations space E. Problem: For every Ac+T, 
and for every g e E, define a subspace Lac such that: (i) Lac is a complement of T in 
E: 

L AC + T = E, L Ac nT = {0}. 



(ii) Ac is the point of minimum of G q on Lac + Ac: 

Ac = arg^nun^ G q (x). (79) 



Besides (i) and (ii), we also impose the requirement of a maximal smoothness (ana- 
lyticity) on Lac as a function of g and Ac. Requirement (79) implies that Ac is the 
quasi-equilibrium point for the given Lac, while the problem in a whole is the inverse 
quasi-equilibrium problem: We construct Lac such that T will be the quasi-equilibrium 
manifold. Then subspaces Lac will actually be the kernels of the quasi-equilibrium pro- 
jector. 

Let f 1 , . . . , f k be the orthonormalized with respect to (•, •) scalar product basis of T, 
vector h be orthogonal to T, (h, h) = 1, g = af 1 + (3h. Condition (79) implies that the 
vector VG is orthogonal to Lac in the point Ac. 

Let us first consider the case (3 = 0. The requirement of analyticity of Lac as the function 
of a and Ac implies Lac — Lq + o(l), where Lq = T 1 - is the orthogonal completement 
of T with respect to scalar product (•,•). The constant solution, Lac = Lq also satisfies 
(79). Let us fix a ^ 0, and extend this latter solution to (3 ^ 0. With this, we obtain a 
basis, h, . . . , l n -k- Here is the simplest construction of this basis: 

/^-(q + AcQfe 

1 (/? 2 + (a + Ac 1 ) 2 ) 1 /2' ^ 



where Aci = (Ac, f ± ) is the first component in the expansion, Ac = J2i ^ c ifi- The rest of 
the basis elements, l 2 , ■ ■ ■ , i n -fc form the orthogonal completement of T© (h) with respect 
to scalar product (•, ■), (h) is the line spanned by h. 

Dependence Lac (80) on Ac, a and (3 is singular: At a + Aci, vector l\ e T, and then 
Lac is not the completement of T in £ anymore. For a^O, dependence Lac gives one 
of the solutions to the inverse quasi-equilibrium problem in the neighborhood of zero in 
T. We are interested only in the limit, 

lim L AC = Lin j^==iJ,Z 2 , ...,l n _ k \ . (81) 



AC 



+ o [V^TJ 2 
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Finally, let us define now the projector Pc of the space E onto T c n. If H^VG G T c fi, 
then P c is the orthogonal projector with respect to the scalar product (•, •): 



k 



Pcz = Y,fi(fn z )- 



(82) 



If H c l VG £ T c £l, then, according to Eq. (81), 



(f 1 ,z)-(l 1 ,z)(f 1 ,l 1 ) 

l-(/l,il> 2 



k 



Pcz 



(83) 



i=2 



where . . . , f k } is the orthonormal with respect to (•, •) basis of T c f2, h is orthogonal 



to T, (h, h) = 1, H c 1 VG = a/ 1 +/3h, Zi = (Pf^aty/y/aTW, </i, *i> = 



Thus, for solving the invariance equation iteratively, one needs only projector P c (83), 
and one does not need a priori parameterization of ft anymore. 



10 Method of invariant grids 

Elsewhere above in this paper, we considered the invariant manifold, and methods for their 
construction, without addressing the question of how to implement it in a constructive way. 
In most of the works (of us and of other people on similar problems), analytic forms were 
needed to represent manifolds. However, in order to construct manifolds of a relatively low 
dimension, grid-based representations of manifolds become a relevant option. The Method 
of invariant grids (MIG) was suggested recently in [32]. 

The principal idea of (MIG) is to find a mapping of finite-dimensional grids into the phase 
space of a dynamic system. That is we construct not just a point approximation of the 
invariant manifold, but an invariant grid. When refined, in the limit it is expected to 
converge, of course, to the invariant manifold, but it is a separate, independently defined 
object. 

Let's denote L = R n , G is a discrete subset of R n . A natural choice would be a regular 
grid, but, this is not crucial from the point of view of the general formalism. For every 
point y G G, a neighborhood of y is defined: V y C G, where V y is a finite set, and, in 
particular, y G V y . On regular grids, V y includes, as a rule, the nearest neighbors of y. It 
may also include next to nearest points. 

For our purposes, one should define a grid differential operator. For every function, defined 
on the grid, also all derivatives are defined: 



dyi yeG 




Z£Vy 



(84) 
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where qi(z,y) are some coefficients. 

Here we do not specify the choice of the functions qi(z, y). We just mention in passing that, 
as a rule, equation (84) is established using some interpolation of / in the neighborhood 
of y in R n by some differentiate functions (for example, polynomial). This interpolation 
is based on the values of / in the points of V y . For regular grids, qi(z,y) are functions of 
the difference z — y. For some ys which are close to the edges of the grid, functions are 
defined only on the part of V y . In this case, the coefficients in (84) should be modified 
appropriately in order to provide an approximation using available values of /. Below we 
will assume this modification is always done. We also assume that the number of points 
in the neighborhood V y is always sufficient to make the approximation possible. This 
assumption restricts the choice of the grids G. Let's call admissible all such subsets G, 
on which one can define differentiation operator in every point. 

Let F be a given mapping of some admissible subset G C R n into U. For every y e V we 
define tangent vectors: 

T y = Linfo}?, (85) 

where vectors g^i = 1, . . .n) are partial derivatives (84) of the vector-function F: 
OF 

9i = -fT = 9i(z,y)F(z), (86) 

^ Z£Vy 



or in the coordinate form: 
OF 

^ = 7)T = E*(*> !/)*>(*)■ ( 87 ) 



Here (g^j is the jih coordinate of the vector (grj, and Fj(z) is the jth coordinate of the 
point F(z). 

The grid G is invariant, if for every node y G G the vector field J(F(y)) belongs to the 
tangent space T y (here J is the right hand site of the kinetic equations (4)). 

So, the definition of the invariant grid includes: 

1) Finite admissible subset G C R n ; 



2) A mapping F of this admissible subset G into U (where U is the phase space for kinetic 
equations (4)); 

3) The differentiation formulas (84) with given coefficients qi(z,y); 
The grid invariance equation has a form of inclusion: 
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J(F(y)) eT y for every y e G, 



or a form of equation: 



(1 - P F{y) )J(F(y)) = for every y e G, 
where Pf( v ) is the thermodynamic projector (83). 

The grid differentiation formulas (84) are needed, in the first place, to establish the tangent 
space T y , and the null space of the thermodynamic projector Pf{ v ) in each node. It is 
important to realise that locality of construction of thermodynamic projector enables this 
without a need for a global parametrization. 

Basically, in our approach, the grid specifics are in: (a) differentiation formulas, (b) grid 
construction strategy (the grid can be extended, contracted, refined, etc.) The invariance 
equations (45), the iteration Newton method (46,47), and the formulas of the relaxation 
approximation (77) do not change at all. For convenience, let us repeat all these formulas 
in the grid context. 

Let c = F(y) be position of a grid's node y immersed into phase space U. We have set 
of tangent vectors g^x), defined in c (86), (87). Thus, the tangent space T y is defined 
by (85). Also, one has the thermodynamic Lyapunov function G(c), the linear functional 
D c G\c, and the subspace T 0y = T y f)ker D c G\c in T y . Let T 0y ^ T y . In this case we have 
a vector e y e T y , orthogonal to T 0y , D c G\c(e y ) = 1- Then, the thermodynamic projector 
is defined as: 

Pc* = Poc • +e y D c G\ c ; (88) 



where P c is the orthogonal projector on T 0y with respect to the entropic scalar product 
(.).- 

If T 0y = T y , then the thermodynamic projector is the orthogonal projector on T y with 
respect to the entropic scalar product (,)c- 

For the Newton method with incomplete linearization, the equations for calculating new 
node position c' = c + 8c are: 

P c Sc = 

(89) 

(1 - P c )(J(c) + DJ(c)6c) = 0. 



Here DJ(c) is a matrix of derivatives of J, calculated in c. The self-adjoint linearization 
may be useful too (see section 5.2). 
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Equation (89) is a system of linear algebraic equations. In practice, it is convenient to 
choose some orthonormal (with respect to the entropic scalar product) basis 6j in kerP c . 
Let r = dim{kerP c ). Then 5c = Yn=i and the system looks like 

r 

J25 k {b i ,DJ(c)b k )c = -(J(c),b i )c,i = l...r. (90) 
k=i 



Here (,) c is the entropic scalar product (28). This is the system of linear equations for 
adjusting the node position accordingly to the Newton method with incomplete lineariza- 
tion. 

For the relaxation method, one needs to calculate the defect A c = (1 — P c ) J(c), and the 
relaxation step 

T ( x ) = (Ac,A c ) c 

T[X) (A c ,DJ(c)A c ) c (91j 



Then, new node position x' is calculated as 

c' = c + r(c)A c . (92) 



This is the equation for adjusting the node position according to the relaxation method. 



10.1 Grid construction strategy 



From all reasonable strategies of the invariant grid construction we will consider here the 
following two: growing lump and invariant flag. 



10.1.1 Growing lump 

In this strategy one chooses as initial the equilibrium point y* . The first approximation is 
constructed as F(y*) = c*, and for some initial Vq (V y * C Vo) one has F(y) = c*+A(y—y*), 
where A is an isometric embedding (in the standard Euclidean metrics) of R n in E. 

For this initial grid one makes a fixed number of iterations of one of the methods chosen 
(Newton's method with incomplete linearization or the relaxation method), and, after 
that, puts V\ = UyeVo Vy an d extends F from V onto V\ using linear extrapolation and 
the process continues. One of the possible variants of this procedure is to extend the grid 
from Vi to Vi+i not after a fixed number of iterations, but when the invariance defect A y 
becomes smaller than a given e (in a given norm, which is entropic, as a rule), for all 
nodes y E Vi. The lump stops growing when it reaches the boundary and is within a given 
accuracy ||A|| < e. 
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10.1.2 Invariant flag 



For the invariant flag one uses sufficiently regular grids G, in which many points are 
situated on the coordinate lines, planes, etc. One considers the standard flag R° C R 1 C 
R 2 C ... C R n (every next space is constructed by adding one more coordinate). It 
corresponds to a succession of grids {y} C G 1 C G 2 ... C G n , where {y*} = R°, and G* 
is a grid in R\ 

First, y* is mapped in c* and further F(y*) = c*. Then an invariant grid is constructed 
on V 1 C G 1 (up to the boundaries U and within a given accuracy ||A|| < e). After the 
neighborhoods in G 2 are added to the points V 1 , and, using such extensions, the grid 
V 2 C G 2 is constructed (up to the boundaries and within a given accuracy) and so on, 
until V n C G" will be constructed. 

We must underline here that, constructing the k-th grid V k C G fc , the important role of 
the grids of smaller dimension V° C ... C V k ^ 1 C V k embedded in it, is preserved. The 
point F(y*) = x* is preserved. For every y e V q (q < k) the tangent vectors gi, ...,g q are 
constructed, using the differentiation operators (84) on the whole V k . Using the tangent 
space T y = Lin{gi, .., g q }, the projector PF{ y ) is constructed, the iterations are applied 
and so on. All this is done to obtain a succession of embedded invariant grids, given by 
the same map F. 

10.1.3 Boundaries check and the entropy 

We construct grid mapping of F onto the finite set V G G. The technique of checking 
if the grid still belongs to the phase space U of kinetic system (F(V) C U) is quite 
straightforward: all the points y <E V are checked to belong to U. If at the next iteration a 
point F(y) leaves U, then it is returned inside by a homothety transform with the center 
in x*. Since the thermodynamic Lyapunov function is a convex function, the homothety 
contraction with the center in x* decreases it monotonously. Another variant is cutting 
off the points leaving U . 

By the way it was constructed (83) the kernel of the thermodynamic projector is annulled 
by the entropy differential. Thus, in the first order, steps in the Newton method with in- 
complete linearization (46,47) as well as in the relaxation methods (76), (77) do not change 
the entropy. But, if the steps are quite large, then the increasing of the thermodynamic 
Lyapunov function can become essential and the points are returned on their level by the 
homothety contraction with the center in the equilibrium point. 

10.2 Instability of fine grids 

When one reduces the grid step (spacing between the nodes) in order to get a finer grid, 
then, starting from a definite step, it is possible to face the problem of the Courant 
instability [16]. Instead of converging, at the every iteration the grid becomes entangled 
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(see Fig. 2). 



The way to get rid off this instability is well-known. This is decreasing the time step. 
Instead of the real time step, we have a shift in the Newtonian direction. Formally we can 
assign for one complete step in the Newtonian direction a value h — 1. Let us consider 
now the Newton method with an arbitrary h. For this, let us find 8c = SF(y) from (89), 
but we will change 5c proportionally to h: the new value of c n+ i = F n+ i(y) will be equal 
to 

F n+1 (y) = F n (y) + h n 5F n (y) (93) 
where the lower index n denotes the step number. 

One way to choose the h step value is to make it adaptive, controlling the average value 
of the invariance defect || A y || at every step. Another way is the convergence control: then 
J2 h n plays a role of time. 

Elimination of Courant instability for the relaxation method can be made quite analo- 
gously. Everywhere the step h is maintained as big as it is possible without convergence 
problems. 

10.3 What space is the most appropriate for the grid construction? 

For the kinetics systems there are two distinguished representations of the phase space: 

• The densities space (concentrations, energy or probability densities, etc.) 

• The spaces of conjugate intensive quantities, potentials (temperature, chemical poten- 
tials, etc.) 

The density space is convenient for the construction of quasi-chemical representations. 
Here the balance relations are linear and the restrictions are in the form of linear inequal- 
ities (the densities themselves or some linear combinations of them must be positive). 

The conjugate variables space is convenient in the sense that the equilibrium conditions, 
given the linear restrictions on the densities, are in the linear form (with respect to the 
conjugate variables). In these spaces the quasi-equilibrium manifolds exist in the form of 
linear subspaces and, vise versa, linear balance equations turns out to be equations of the 
conditional entropy maximum. 

The duality we've just described is very well-known and studied in details in many works 
on thermodynamics and Legendre transformations [13,73]. In the previous section, the 
grids were constructed in the density space. But the procedure of constructing them in 
the space of the conjugate variables seems to be more consistent. The principal argu- 
ment for this is the specific role of quasi-equilibrium, which exists as a linear manifold. 
Therefore, linear extrapolation gives a thermodynamically justified quasi-equilibrium ap- 
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proximation. Linear approximation of the slow invariant manifold in the neighborhood of 
the equilibrium in the conjugate variables space already gives the global quasi-equilibrium 
manifold, which corresponds to the motion separation (for slow and fast motions) in the 
neighborhood of the equilibrium point. 

For the mass action law, transition to the conjugate variables is simply the logarithmic 
transformation of the coordinates. 



10. 4 Carleman's formulas in the analytical invariant manifolds approximations. First 
profit from analyticity: super-resolution 

When constructing invariant grids, one must define the differential operators (84) for 
every grid node. For calculating the differential operators in some point y, an interpolation 
procedure in the neighborhood of y is used. As a rule, it is an interpolation by a low- 
order polynomial, which is constructed using the function values in the nodes belonging 
to the neighbourhood of y in G. This approximation (using values in the closest nodes) 
is natural for smooth functions. But, for the systems (4) with analytical right hand side 
we are looking for the analytical invariant manifold (due to Lyapunov auxiliary theorem 
[60,54]). Analytical functions have much more "rigid" structure than the smooth ones. 
One can change a smooth function in the neighborhood of any point in such a way, that 
outside this neighborhood the function will not change. In general, this is not possible for 
analytical functions: a kind of "long-range" effect takes place (as is well known) . 

The idea is to use this effect and to reconstruct some analytical function /g using function 
given on G. There is one important requirement: if these values on G are values (given in 
the points of G) of some function / which is analytical in the given neighborhood U, then 
if the G is refined "correctly", one must have f G — > /. The sequence of reconstructed 
function / G should should converge to the "proper" function /. 

What is the "correct refinement"? For smooth functions for the convergence /g — > / it 
is necessary and sufficient that, in the course of refinement, G would approximate the 
whole U with arbitrary accuracy. For analytical functions it is only necessary that, under 
the refinement, G would approximate some uniqueness set 3 A C U . Suppose we have 
a sequence of grids G, each next is finer than previous, which approximates a set A. 
For smooth functions, using function values defined on the grids, one can reconstruct the 
function in A. For analytical functions, if the analyticity area U is known, and A is a 
uniqueness set in U, then one can reconstruct the function in U . The set U can be essen- 
tially bigger than A; because of this such extension was named as superresolution effects 
[1,38]. There exist constructive formulas for construction of analytical functions / G for 
different areas U, uniqueness sets A C U and for different ways of discrete approximation 
of A by a sequence of fined grids G [1]. Here we provide only one Carleman's formula 
which is the most appropriate for our purposes. 

3 Let's remind to the reader that A C U is called uniqueness set in U if for analytical in U 
functions ip and (p from iJj\a = ¥>\a it follows ip = ip. 
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Let area U — Q™ C C n be a product of strips Q a C C, Q a = {z\lmz < a}. We will 
construct functions holomorphic in Q™. This is effectively equivalent to the construction 
of real analytical functions / in whole R n with a condition on the convergence radius 
r(x) of the Taylor series for / as a function of each coordinate: r(x) > a in every point 
x G R n . 

The sequence of fined grids is constructed as follows: let for every I = 1, ...,n a finite 
sequence of distinct points N t C D a be defined: 

N i = i x ij\j = 1,2,3...},^- ^ xu for i ^ j (94) 

The uniqueness set A, which is approximated by a sequence of fined finite grids, has the 
form: 

A = JVi x N 2 x ... x N n = {(x lil ,x 2 i 2 ,..,x nin )\i 1 ^ n = 1,2,3,...} (95) 

The grid G m is defined as the product of initial fragments Ni of length m: 

G m = {(xi il ,x 2 i 2 ...x nin )\l < ii v . >n < m} (96) 



Let's denote A = 2a /n (a is a half- width of the strip Q a ). The key role in the construction 
of the Carleman's formula is played by the functional cj^(-u,p, Z) of 3 variables: u e U = 
Q™, p is an integer, 1 < p < m, I is an integer, 1 < p < n. Further u will be the coordinate 
value in the point where the extrapolation is calculated, / will be the coordinate number, 
and p will be an element of multi-index {ii, ...,i n } for the point (xu 1 ,X2i 2 , ■ ■■,x n i n ) £ G: 

u m (u,p,l) - - e x Slp ^ u _ Xlp y*xi v X jj (e A ^ - e AiE y)(e Au + e Ai y) ^ 

For real-valued formula (97) becomes simpler: 

cu A ( M ,p,/) = 2 w - , 6 W -x TT r> > x ( 98 ) 

The Carleman's formula for extrapolation from Gm on U — Q% (a = 7rA/2) has the form 
(z = (z 1 ,...,z n )): 

m n 

Uz)= £ fMU^i^jJ), (99) 

where k = k u ..,k n , x k = (x lkl , x 2k2 , x nkn ). 
There exists a theorem [1]: 
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If / e H 2 (Q™), then f(z) = lim m ^ocf m (z), where H 2 (Q™) is the Hardy class of 
holomorphic in Q™ functions. 

It is useful to present the asymptotics of (99) for big |Rez,-|. For this we will consider the 
asymptotics of (99) for big |Reu|: 



2 m e *X lp _|_ gAXy 

n 



+ o(|Rew|- 1 ). 



(100) 



From the formula (99) one can see that for the finite m and |ReZj| — > oo function \ f m {z)\ 
behaves like const • f\j 

This property (null asymptotics) must be taken into account when using the formula 
(99). When constructing invariant manifolds F(W), it is natural to use (99) not for the 
immersion F(y), but for the deviation of F(y) from some analytical ansatz Fo(y). 

The analytical ansatz F (y) can be obtained using Taylor series, just as in the Lyapunov 
auxiliary theorem [60]. Another variant is using Taylor series for the construction of Pade- 
approximations. 

It is natural to use approximations (99) in dual variables as well, since there exists for 
them (as the examples demonstrate) a simple and very effective linear ansatz for the 
invariant manifold. This is the slow invariant subspace E s \ ow of the operator of linearized 
system (4) in dual variables in the equilibrium point. This invariant subspace corresponds 
to the the set of "slow" eigenvalues (with small |ReA|, ReA < 0). In the initial space (of 
concentrations or densities) this invariant subspace is the quasi-equilibrium manifold. It 
consist of the maximal entropy points on the afline manifolds of the x + £fast form, where 
-Efast is the "fast" invariant subspace of the operator of linearized system (4) in the initial 
variables in the equilibrium point. It corresponds to the "fast" eigenvalues (big |ReA|, 
ReA < 0). 



10.5 Example: Two-step catalytic reaction 

Let's consider a two-step four- component reaction with one catalyst A 2 : 

A 1 + A 2 <-> A 3 <-> A 2 + A 4 (101) 



We assume the Lyapunov function of the form G = J2t=i Ci[ln(ci/c*) — 1]. The kinetic 
equation for the four-component vector of concentrations, c = (ci, c 2 , C3, C4), has the form 

c = ll W 1 + l2 W 2 . (102) 
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Here 71,2 are stoichiometric vectors, 

7l = (-1,-1, 1,0), 72 = (0,1, -1,1), 



(103) 



while functions W^ 2 are reaction rates: 



Wi = k\c x c 2 - k l c 3 , W 2 = k£c 3 - k 2 c 2 c A . 



(104) 



Here k 12 are reaction rate constants. The system under consideration has two conservation 
laws, 



or (61,2, c) = B 12 , where 61 = (1, 0, 1, 1) and bi = (0, 1, 1, 0). The nonlinear system (101) 
is effectively two-dimensional, and we consider a one-dimensional reduced description. For 
our example, we chose the following set of parameters: 



B 1 = 1.0, B 2 = 0.2 

In Fig. 3 one-dimensional invariant grid is shown in the (ci,c 4 ,c 3 ) coordinates. The grid 
was constructed by growing the grid, as described above. We used Newtonian iterations 
to adjust the nodes. The grid was grown up to the boundaries of the phase space. 

The grid derivatives for calculating tangent vectors g were taken as simple as g(ci) = 
(cj + i-Cj_i)/||cj+i-Cj_i|| for the internal nodes and gr(ci) = (c 1 -c 2 )/\\c 1 -c 2 \\, g(Cn) = 
(c„ — Cn-i)/\\c n — c„_i|| for the grid's boundaries. Here Xj denotes the vector of the ith 
node position, n is the number of nodes in the grid. 

Close to the phase space boundaries we had to apply an adaptive algorithm for choosing 
the time step h: if, after the next growing step and applying N = 20 complete Newtonian 
steps, the grid did not converged, then we choose a new h n+ i = h n /2 and recalculate the 
grid. The final value for h was h ~ 0.001. 

The nodes positions are parametrized with entropic distance to the equilibrium point 
measured in the quadratic metrics given by He = \\d 2 G(c) /dcidcjW in the equilibrium 
c*. It means that every node is on a sphere in this quadratic metrics with a given radius, 
which increases linearly. On this figure the step of the increase is chosen to be 0.05. 
Thus, the first node is on the distance 0.05 from the equilibrium, the second is on the 
distance 0.10 and so on. Fig. 4 shows several basic values which facilitate understanding 
of the object (invariant grid) extracted. The sign on the x-axis of the graphs at Fig. 4 
is meaningless, since the distance is always positive, but in this situation it denotes two 
possible directions from the equilibrium point. 



ci + c 3 + c 4 = B 1 , c 2 + c 3 = B- 



2, 



(105) 



kj = 0.3, ki = 0.15, kj = 0.8, k 2 = 2.0; 
c * = 0.5, c* 2 = 0.1, c* 3 = 0.1, c\ = 0.4; 



(106) 
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Fig. 4a,b effectively represents the slow one- dimensional component of the dynamics of 
the system. Given any initial condition, the system quickly finds the corresponding point 
on the manifold and starting from this point the dynamics is given by a part of the graph 
on the Fig. 4a,b. 

One of the useful values is shown on the Fig. 4c. It is the relation between the relaxation 
times "toward" and "along" the grid (A2/A1, where Ai,A2 are the smallest and the second 
smallest by absolute value non-zero eigenvalue of the system, symmetrically linearized in 
the point of the grid node). It shows that the system is very stiff close to the equilib- 
rium point, and less stiff (by one order of magnitude) on the borders. This leads to the 
conclusion that the reduced model is more adequate in the neighborhood of the equi- 
librium where fast and slow motions are separated by two orders of magnitude. On the 
very end of the grid which corresponds to the positive absciss values, our one-dimensional 
consideration faces with definite problems (slow manifold is not well-defined). 

10.6 Example: Model hydrogen burning reaction 

In this section we consider a more interesting illustration, where the phase space is 6- 
dimensional, and the system is 4-dimensional. We construct an invariant flag which con- 
sists of 1- and 2-dimensional invariant manifolds. 

We consider chemical system with six species called (provisionally) H 2 (hydrogen), 2 
(oxygen), H 2 (water), H, O, OH (radicals). We assume the Lyapunov function of the 
form G = Yh=i c S n ( c i/ c i) — !]• The subset of the hydrogen burning reaction and corre- 
sponding (direct) rate constants have been taken as: 



1. 


H 2 2H 


kf = 


2 


2. 


2 «-> 20 


k 2 — 


1 


3. 


H 2 <-> H + OH 


kf = 


1 


4. 


H 2 + <-> H + OH 


&4~ = 


10 3 


5. 


2 + H ^O + OH 


= 


10 3 


6. 


H 2 + <-> H 2 




10 2 



(107) 



The conservation laws are: 

2ch 2 + 2c H , 2 o + c H + c H = b H 

(108) 

2c 02 + c H20 + c + c Q H = b 

For parameter values we took bu = 2, bo = 1, and the equilibrium point: 

c* = 0.27 c* 09 = 0.135 c* mo = 0.7 c* H = 0.05 c* = 0.02 c* OH = 0.01 (109) 
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Other rate constants k { ,i — 1..6 were calculated from c* value and kf . For this system 
the stoichiometric vectors are: 

7i = (-1, 0, 0, 2, 0, 0) 72 = (0, -1, 0, 0, 2, 0) 

73 = (0,0, -1,1, 0,1) 7 4 = (-1,0, 0,1, -1,1) (110) 
7 5 = (0, -1, 0, -1, 1, 1) 7e = (-1, 0, 1, 0, -1, 0) 



We stress here once again that the system under consideration is fictional in that sense 
that the subset of equations corresponds to the simplified picture of this physical-chemical 
process and the constants do not correspond to any measured ones, but reflect only basic 
orders of magnitudes of the real-world system. In this sense we consider here a qualitative 
model system, which allows us to illustrate the invariant grids method without excessive 
complication. Nevertheless, modeling of real systems differs only in the number of species 
and equations. This leads, of course, to computationally harder problems, but not the 
crucial ones, and the efforts on the modeling of real-world systems are on the way. 

Fig. 5a presents a one-dimensional invariant grid constructed for the system. Fig. 5b shows 
the picture of reduced dynamics along the manifold (for the explanation of the meaning of 
the x-coordinate, see the previous subsection). On Fig. 5c the three smallest by absolute 
value non-zero eigen values of the symmetrically linearized system A sym have been shown. 
One can see that the two smallest values "exchange" on one of the grid end. It means 
that one-dimensional "slow" manifold has definite problems in this region, it is just not 
defined there. In practice, it means that one has to use at least two-dimensional grids 
there. 



Fig. 6a gives a view onto the two-dimensional invariant grid, constructed for the system, 
using the "invariant flag" strategy. The grid was grown starting from the ID-grid con- 
structed at the previous step. At the first iteration for every node of the initial grid, two 
nodes (and two edges) were added. The direction of the step was chosen as the direction 
of the eigenvector of the matrix A sym (in the point of the node), corresponding to the 
second "slowest" direction. The value of the step was chosen to be e = 0.05 in terms of 
entropic distance. After several Newtonian iterations done until convergence, new nodes 
were added in the direction "ortogonal" to the ID-grid. This time it is done by linear 
extrapolation of the grid on the same step e = 0.05. When some new nodes have one or 
several negative coordinates (the grid reaches the boundaries) they were cut off. If a new 
node has only one edge, connecting it to the grid, it was excluded (since it does not allow 
calculating 2D-tangent space for this node). The process continues until the expansion is 
possible (after this, every new node has to be cut off). 

Strategy of calculating tangent vectors for this regular rectangular 2D-grid was chosen to 
be quite simple. The grid consists of rows, which are co-oriented by construction to the 
initial ID-grid, and columns that consist of the adjacent nodes in the neighboring rows. 
The direction of "columns" corresponds to the second slowest direction along the grid. 
Then, every row and column is considered as ID-grid, and the corresponding tangent 
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vectors are calculated as it was described before: 




(Cfc,i+l — Ck,i-l) /\\Ck,i+l — Ck,i-l 



for the internal nodes and 



9ro W ( c k,l) = (Cfc,l - C fc)2 )/||c fc> i 



c fc,n fe -l)/ || c fc,rtfc Cfc infc _i 



for the nodes which are close to the grid's edges. Here Xk,i denotes the vector of the node in 
the kth row, ith column; is the number of nodes in the fcth row. Second tangent vector 
9coi( c k,i) is calculated completely analogously. In practice, it is convenient to orthogonalize 
9row( c k,i) and g co i{c kA ). 

Since the phase space is four-dimensional, it is impossible to visualize the grid in one of the 
coordinate 3D-views, as it was done in the previous subsection. To facilitate visualization 
one can utilize traditional methods of multi-dimensional data visualization. Here we make 
use of the principal components analysis (see, for example, [45]), which constructs a three- 
dimensional linear subspace with maximal dispersion of the othogonally projected data 
(grid nodes in our case). In other words, method of principal components constructs in 
multi-dimensional space such a three-dimensional box inside which the grid can be placed 
maximally tightly (in the mean square distance meaning). After projection of the grid 
nodes into this space, we get more or less adequate representation of the two-dimensional 
grid embedded into the six-dimensional concentrations space (Fig. 6b). The disadvantage 
of the approach is that the axes now do not have explicit meaning, being some linear 
combinations of the concentrations. 

One attractive feature of two-dimensional grids is the possibility to use them as a screen, 
on which one can display different functions /(c) defined in the concentrations space. This 
technology was exploited widely in the non-linear data analysis by the elastic maps method 
[41]. The idea is to "unfold" the grid on a plane (to present it in the two-dimensional 
space, where the nodes form a regular lattice). In other words, we are going to work in 
the internal coordinates of the grid. In our case, the first internal coordinate (let's call 
it si) corresponds to the direction, co-oriented with the one- dimensional invariant grid, 
the second one (let's call it s 2 ) corresponds to the second slow direction. By how it was 
constructed, s 2 = line corresponds to the one-dimensional invariant grid. Units of s± 
and S2 are entropic distances in our case. 

Every grid node has two internal coordinates (s 1 , s 2 ) and, simultaneously, corresponds to a 
vector in the concentration space. This allows us to map any function /(c) from the multi- 
dimensional concentration space to the two-dimensional space of the grid. This mapping 
is defined in a finite number of points (grid nodes), and can be interpolated (linearly, 
in the simplest case) in between them. Using coloring and isolines one can visualize the 
values of the function in the neighborhood of the invariant manifold. This is meaningful, 
since, by the definition, the system spends most of the time in the vicinity of the invariant 
manifold, thus, one can visualize the behaviour of the system. As a result of applying the 
technology, one obtains a set of color illustrations (a stack of information layers), put 
onto the grid as a map. This allows applying all the methods, working with stack of 
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information layers, like geographical information systems (GIS) methods, which are very 
well developed. 

In short words, the technique is a useful tool for exploration of dynamical systems. It 
allows to see simultaneously many different scenarios of the system behaviour, together 
with different system's characteristics. 

The simplest functions to visualize are the coordinates: q(c) = q. On Fig. 7 we displayed 
four colorings, corresponding to the four arbitrarily chosen concentrations functions (of 
H 2 , O, H and OH; Fig. 7a-d). The qualitative conclusions that can be made from the 
graphs are that, for example, the concentration of H 2 practically does not change during 
the first fast motion (towards the ID-grid) and then, gradually changes to the equilibrium 
value (the H 2 coordinate is "slow"). The O coordinate is the opposite case, it is "fast" 
coordinate which changes quickly (on the first stage of motion) to the almost equilibrium 
value, and then it almost does not change. Basically, the slope angles of the coordinate 
isolines give some presentation of how "slow" a given concentration is. Fig. 7c shows 
interesting behaviour of the OH concentration. Close to the ID grid it behaves like "slow 
coordinate" , but there is a region on the map where it has clear "fast" behaviour (middle 
bottom of the graph). 

The next two functions which one can want to visualize are the entropy S = —G and the 
entropy production <r(c) = —dG/dt(c) = — J2i hi(q/c*)Q. They are shown on Fig. 8a,b. 

Finally, we visualize the relation between the relaxation times of the fast motion towards 
the 2D-grid and along it. This is given on the Fig. 8c. This picture allows to make a con- 
clusion that two-dimensional consideration can be appropriate for the system (especially 
in the "high H 2 , high O" region), since the relaxation times "towards" and "along" the 
grid are definitely separated. One can compare this to the Fig. 8d, where the relation 
between relaxation times towards and along the ID-grid is shown. 



11 Method of invariant manifold for open systems 



One of the problems to be focused on when studying closed systems is to prepare exten- 
sions of the result for open or driven by flows systems. External flows are usually taken 
into account by additional terms in the kinetic equations (4): 

c = j(c) + n. (111) 



Zero-order approximation assumes that the flow does not change the invariant manifold. 
Equations of the reduced dynamics, however, do change: Instead of J(c(M)) we substitute 
J(c(M)) + n into Eq. (39): 

Mi = (VM,| C(M) , J(c(M)) + n). (112) 
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Zero-order approximation assumes that the fast dynamics in the closed system strongly 
couples the variables c, so that flows cannot influence this coupling. 

First-order approximation takes into account the shift of the invariant manifold by 5c. 
Equations for Newton's iterations have the same form (57) but instead of the vector field 
J they take into account the presence of the flow: 



where projector Pc corresponds to the unperturbed manifold. 

The first-order approximation means that fluxes change the coupling between the variables 
(concentrations). It is assumed that these new coupling is also set instantaneously (neglect 
of inertia). 

Remark. Various realizations of the first-order approximation in physical and chemical dy- 
namics implement the viewpoint of an infinitely small chemical reactor driven by the flow. 
In other words, this approximation is applicable in the Lagrangian system of coordinates 
[51,79]. Transition to Eulerian coordinates is possible but the relations between concen- 
trations and the flow will change its form. In a contrast, the more simplistic zero-order 
approximation is equally applicable in both the coordinate system, if it is valid. 



12 Conclusion 

In this paper, we have presented the method for constructing the invariant manifolds 
for reducing systems of chemical kinetics. Our approach to computations of invariant 
manifolds of dissipative systems is close in spirit to the Kolmogorov-Arnold-Moser theory 
of invariant tori of Hamiltonian systems [5,6]: We also base our consideration on the 
Newton method instead of Taylor series expansions [7], and systematically use duality 
structures. Recently, a version of an approach based on the invariance equations was 
rediscovered in [54]. He was solving the invariance equation by a Taylor series expansion. 
A counterpart of Taylor series expansions for constructing the slow invariant manifolds in 
the classical kinetic theory is the famous Chapman-Enskog method. The question of how 
this compares to iteration methods was studied extensively for certain classes of Grad 
moment equations [31,50,48]. 

The thermodynamic parameterization and the selfadjoint linearization arise in a natural 
way in the problem of finding slowest invariant manifolds for closed systems. This also 
leads to various applications in different approaches to reducing the description, in partic- 
ular, to a thermodynamically consistent version of the intrinsic low-dimensional manifold, 
and to model kinetic equations for lifting the reduced dynamics. Use of the thermody- 
namic projector makes it unnecessary global parameterizations of manifolds, and thus 
leads to computationally promising grid-based realizations. 



[1 - P C ](II + L' c 5c) = 0, P c 5c = 0, 
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Invariant manifolds are constructed for closed space-independent chemical systems. We 
also describe how to use these manifolds for modeling open and distributed systems. 
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Fig. 1. Images of the initial quasi-equilibrium manifold (bold line) and the first two corrections 
(solid normal lines) in the phase plane [01,03] for two-step catalytic reaction. Dashed lines are 
individual trajectories. 

Fig. 2. Grid instability. For small grid steps approximations in the calculation of grid derivatives 
lead to the grid instability effect. On the figure several successive iterations of the algorithm 
without adaptation of the time step are shown that lead to undesirable "oscillations", which 
eventually destruct the grid starting from one of it's ends. 

Fig. 3. One-dimensional invariant grid (circles) for two-dimensional chemical system. Projection 
into the 3d-space of ci, C4, C3 concentrations. The trajectories of the system in the phase space 
are shown by lines. The equilibrium point is marked by square. The system quickly reaches the 
grid and further moves along it. 

Fig. 4. One-dimensional invariant grid for two-dimensional chemical system, a) Values of the con- 
centrations along the grid, b) Values of the entropy (— G) and the entropy production {—dG/dt) 
along the grid, c) Relation of the relaxation times "toward" and "along" the manifold. The nodes 
positions are parametrized with entropic distance measured in the quadratic metrics given by 
H c = \\d 2 G(c)/dcidcj\\ in the equilibrium c*. Zero corresponds to the equilibrium. 

Fig. 5. One-dimensional invariant grid for model hydrogen burning system, a) Projection into 
the 3d-space of ch, co, coh concentrations, b) Concentration values along the grid, c) three 
smallest by absolute value non-zero eigen values of the symmetrically linearized system. 

Fig. 6. Two-dimensional invariant grid for the model hydrogen burning system, a) Projection 
into the 3d-space of ch, co, coh concentrations, b) Projection into the principal 3D-subspace. 
Trajectories of the system are shown coming out from the every grid node. Bold line denotes 
the one-dimensional invariant grid, starting from which the 2D-grid was constructed. 

a) Concentration H2 b) Concentration O 

c) Concentration OH d) Concentration H 

Fig. 7. Two-dimensional invariant grid as a screen for visualizing different functions defined 
in the concentrations space. The coordinate axes are entropic distances (see the text for the 
explanations) along the first and the second slowest directions on the grid. The corresponding 
ID invariant grid is denoted by bold line, the equilibrium is denoted by square. 

a) Entropy b) Entropy Production 

c) A3/A2 relation d) A2/A1 relation 

Fig. 8. Two-dimensional invariant grid as a screen for visualizing different functions defined 
in the concentrations space. The coordinate axes are entropic distances (see the text for the 
explanations) along the first and the second slowest directions on the grid. The corresponding 
ID invariant grid is denoted by bold line, the equilibrium is denoted by square. 



53 



This figure "figl.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig2.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig3.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig4.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig5.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig6.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig7.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



This figure "fig8.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/cond-mat/0307076vl 



